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SUMMARY 


A theoretical and experimental investigation has been made of 
the interaction of waves of finite amplitude in ducts of varying 
cross section. A new one-dimensional method is described by 
which the steady-state conditions after passage of an isentropic 
wave may be obtained from a chart without necessity of previ- 
ously required iterative procedures. The chart has proved useful 
in the discussion of the spectrum of physical solutions which has 
been presented herein 

Experiments were carried out in a shock tube to indicate the 
interaction processes 


actual two-dimensional 


to determine how accurately one-dimensional 


nature of the 
occurring, and 
methods could be applied to such interactions. A two-dimen- 
sional channel of area ratio 0.504 was made in such a way that the 
normal shock wave could be incident on either end. Density 
fields were measured with the aid of a Mach-Zehnder inter 
ferometer 

Comparison of the transient density distribution obtained in 
the center of the channel with results of a one-dimensional un 
steady flow calculation for the converging channel indicates ex- 
cellent agreement. The establishment of steady flow was ob- 
served for a wide variety of incident shock strengths using both 
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the converging and the diverging models. The isentropic theory 
weak shocks 


agreement is obtained for strong shocks using real shock theory 


fits the measured densities well for Satisfactory 


INTRODUCTION 


HE PHENOMENA ASSOCIATED with propagation of 
"Was of finite amplitude in channels of varying 
area are of practical as well as fundamental interest in- 
asmuch as such interactions occur in piston engines, 
intermittent jet engines, unsteady flow periods in the 
so-called ‘‘steady flow’ jet engines (for example, pulsat 
ing combustion and supersonic diffuser oscillations in the 
ram-jet), in the starting of supersonic wind tunnels, and 
so on. 

The propagation of a wave of finite amplitude through 
a channel of slowly varying cross section has been 
treated by Guderley! and others.” * The problem can 
be solved by using the one-dimensional method of 
characteristics; the solution is obtained with a graphical 
or numerical calculation. The results of such a calcula- 
tion will be the transient values of the flow quantities. 
The ultimate steady values of the flow properties can be 
obtained by representing the variable area channel by 
a discontinuous step in area and calculating this simple 
wave interaction problem. Schultz-Grunow‘ has pre- 
sented a graphical technique for the calculation of this 
interaction for isentropic flow and waves, whereas a 
numerical method has been given by Rudinger and 
Rinaldi’ for both isentropic and real shock waves. All 
the above methods are of the iterative type. 

An experimental shock tube study of the propagation 
of a shock wave through a channel of varying area has 
been reported by Herzberg and Kantrowitz.6 The 
space-time history of the initial shock in its passage 
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through the channel was compared to the results of a 
one-dimensional characteristics calculation, with good 
agreement being obtained. Bull’ has investigated the 
starting processes in a blowdown wind tunnel when the 
flow is initiated by puncturing a cellophane diaphragm. 

In the present paper, a new method is given for ob- 
taining the steady-state conditions after wave-variable 
area channel interaction for isentropic flow and waves. 
The method features a chart which not only eliminates 
the iterative procedures of the previous methods, but 
which also assists in the discussion of the spectrum of 
A discussion of the nature of the solutions 
The results 


solutions. 
for various initial conditions is presented. 
of an experimental shock tube investigation of the inter- 
action of a shock wave with variable area channels are 
also presented. The experiments consisted of inter- 
ferometric studies of the propagation of shocks of 
varying strengths through a converging and diverging 
channel of 0.504 area ratio. The time required for the 
ultimate steady flow to be established was determined 
in the converging channel for incident shock pressure 
ratios of 1.96 and 3.95. 

The steady-state densities were obtained for both 
channels through a range of incident shock pressure 
ratios; these results are compared to theoretically ob- 
tained values. The transient densities for the con- 
verging channel shock pressure ratio of 1.96 series of 
tests are compared to the results of one-dimensional 
characteristics method calculations. 


SYMBOLS 


A = area 
a = velocity of sound 
L = length of area-change section; 2.85 in. 
M = Mach Number 
p = pressure 
F = Riemann variable; 2a/(y — 1) + u 
Q = Riemann variable; 2a/(y — 1) — u 
R = gas constant 
r = temperature; °R. 
t = time 
u = particle velocity 
x = distance along channel 
y = distance from wall 
Y = height of channel 
Ws = wave moving downstream with respect to fluid 
W- = wave moving upstream with respect to fluid 
7 = ratio of specific heats 
6 = fringe shift 
boundary layer thickness 
p = density 


THEORY 


The ultimate steady-state conditions for wave inter- 
action with a channel of variable area may be simply 
obtained by consideration of the interaction of a wave 
with a channel of discontinuous area change. The 
transient interaction of a wave with a channel of vari- 
able area has also been calculated using this method‘ by 
simulating the channel by a series of constant area tubes 
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connected by discontinuous area steps. A method js 
presented herein for calculating the interaction of 
wave and area discontinuity (for isentropic flow and 
waves) which does not require the iterative procedures 
of the previous methods.* ® A discussion of the nature 
of the solutions for a range of initial and boundary con 


ditions is also presented. 


Wave-Area Change Interaction Chart 


A simple method for obtaining the solution for the 
interaction of waves and area discontinuities as shown 
in Fig. 1 will be developed. Although usually one in- 
cident wave is encountered in practice, the more general 
case of two oppositely propagating incident waves 
striking the area change at the same time will be con- 
sidered. The flow is assumed isentropic and one dimen- 
sional (dependent on distance coordinate x and time { 
only) and all waves are assumed isentropic. For 
simplicity, bands of compression waves or expansion 
waves will be simulated in most of the figures by single 
mean waves. It will be assumed that the reflections of 
the incident waves consist of a single wave in each con- 
stant area section. (Other types of reflections for which 
the chart developed does not give solutions, will be 
pointed out in the next section.) The positive x diree- 
tion will be called the downstream direction. Waves 
which propagate with the stream will be called down- 
stream propagating waves and will be designated |W. 
Waves which propagate against the stream will be called 
upstream propagating waves and will be designated 
W. 

Consider now the interaction indicated in Fig. 1 in 
which a duct of variable area is simulated by an area 
discontinuity. The velocity and speed of sound at sta- 
tions 3, 4, 5, and 6 and the area ratio A,/A» are known; 
flow conditions at stations 1 and 2 are to be deter- 
mined. Since the flow from | to 2 is steady and isen- 
tropic, these states are related by the following equa- 


tions: 
Energy equation: 
uy? + [2/(y — 1) Jay? = wo? + [2/(y — 1) Ja? (J 
Continuity equation: 
pitiA, = prltrA, 


Using the isentropic relation p ~ p’, the equation of 
state p = pRT, and the equation for the speed of sound 


a? = 7RT 
then 


p2/ pr = (de a)" 


which combined with the continuity equation, Eq. (2), 


yields 
1 l 2/(4 1 
‘ A, = Uvlly ‘ A» {oO 


9 
4,Q, 
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ScHEME OF Wave INTERACTIONS AT AW 


AREA DiscOonTINUITY (WV A CHANNEL 


States | and 3 and states 2 and 5 can be related by 
using the Riemann invariant equations.* The Riemann 
equations that apply to finite amplitude isentropic 
waves in a constant area tube, state that the invariants 
Pand Q are constant across upstream (IV’_) and down- 
stream (I’,) propagating waves, respectively, where 


—l)ja+u 


P= [2/(y¥ 
= [2/(y — l)la — u 


) 


a. 


Other authors have used the symbols 7, s, and A, u for 


P. Q.) 


Thus 


P, = Px: [2/(y — 1) Jay + m = [2/(y — 1) Jaz; +4 
(4) 


v2 J) [2 ‘= 1) Jas — w= 
[2 yas 1) Jas — (5) 


Eqs. (1), (3), (4), and (5) constitute a set of four alge- 
braic equations in the four unknowns 1%, d), Ms, and dp. 
Elimination of any three unknown yields, however, 
a polynomial of higher order than three, which, there- 
fore, cannot be solved analytically.* 

* The graphical iterative method of Schultz-Grunow utilizes 
the w — a plane on which are plotted the energy ellipses, Eq. (1): 
— 1)]a? = constant 


u? + [2/(4 


T vA 


and mass flow per unit area hyperbolas, Eq. (3): 


19 + 1 
eae? = constant 


The straight lines representing Eqs. (4) and (5) are drawn through 


the known points (a3, a3) and (#5, a5). The desired solution must 


The combination of Eqs. (1) and (3) yields 


A. M, §1 + [(y — 1)/2]My"" 
A, Me \1 + [(y — 1)/2]Mm9 
where the Mach Number .J/ equals w/a. Division of 

Eq. (4) by Eq. (5) gives 
P, P, j/; +- [2 is = 1) ay 
VI: WN as 


From Eq. (1) it can be determined that 


(“) 
de 
Introduction of the above relation into Eq. (7 
P, Ps [2/(y — lI] + M 
OQ OO; L2/(y¥-)] - 
Fe + [2/(y -— 1 | 
\ 

M,? aT {2 & ieee 1) | 

It will be noted that Eq. (6) gives the boundary con- 
ditions A»/A,f as a function of 1/7, and J/:, whereas Eq. 


M+ [2/(y — 1)] 
M,? + [2/(7 -— 1)] 


gives 


lie on these straight lines |i.e., (1, a) lies on the line through 
(uz, a3), ete.], the states 1 and 3 lying on the same energy ellipse, 
and on two hyperbolas which satisfy the continuity, Eq. (3 

+ A second boundary condition may exist which may be of im 
portance in determining a solution. The duct connecting sections 
A, and A» may have a nonmonotonic variation in the cross-section 
area. Thus a throat may be present, which may or may not de 


mand a sonic throat velocity 
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I. 
J 
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M, 
r=14 
Fic. 2. Chart for solution of wave interaction problems at an 
area discontinuity in a channel 
(8) gives the initial conditions P;/Q; = P/Q: as a 


function of AJ; and A/,. Thus, both Eqs. (6) and (8) 
can be plotted on the J4/,, A/, plane yielding curves of 
constant A»,/A,; and P;/Q2 = P3/Qs5 (see Fig. 2). The 
curves need only be plotted for J/; and A/, > 0, since 
M, and M, must have the same sign. For given area 
ratio and initial conditions P,/Q. = P:;/Q; the final 
steady-state Mach Numbers ./,; and ./, may be deter- 
mined by the intersections of the A2,/A,; and P\/Q: = 
P;/Qs curves. Thus, once a set of curves as illustrated 
in Fig. 2 is constructed, the solution to the interaction 
problem can be obtained quite simply.* 

It will be noted that the chart gives /, as a function 
of M,, with As/A; and P,/Q», as parameters, for all 
The curve OG (A,/A, = 1) 
The curve AI 


isentropic steady flows.* 
gives the trivial steady flow J; = M2. 
(A»/A; = 1) gives the isentropic stationary compres- 
The curve EA (A2/A; = 1) gives what might be 
the expansion” solu- 


sion. 
called 
tion. 


“isentropic stationary 


Discussion of Solutions 


A glance at the interaction chart shown in Fig. 2 re- 
veals that in general the lines of constant P;/Q» cross the 
constant A»/A, lines twice, thus indicating two solu- 
tions for the interaction. In addition it will be noted 
that in certain regions some P/Q» curves do not inter- 
sect with some A,/A, curves, i.e., no solution is given by 
the chart. A discussion of which of the two solutions is 
the physical solution and of solutions which are not 


* Note that for reflections other than that assumed (i.e., a 
W~ wave in the A; section and a Ws wave in the A» section) 


P,/Q2 is in general not equal to P;/Q3. 2 
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given by the chart will be presented in the following 
paragraphs. 

Similar questions of uniqueness of solutions have 
arisen in gas dynamical problems, a notable one arising 
in oblique shock theory. Such uniqueness questions 
are discussed by Courant and Friedrichs;* it appears 
that rigorous mathematical uniqueness arguments of 
these problems are in a primitive state. To test the 
solution obtained for any given example, one could 
calculate the equilibrium solution using the unsteady 
one-dimensional method of characteristics. To insure 
that the calculated solution is the physically unique 
solution, one would have to determine whether the 
solution is stable. This could be done by sending smal] 
amplitude wave pulses through the channel and study- 
ing its effect on the established flow. It appears that if 
the equilibrium flow is entirely subsonic, such a small 
disturbance could not cause a breakdown of the flow, 
The results of Kantrowitz’ and Hess!’ could be used for 
stability studies of equilibrium flows with supersonic 
régimes. 

Inasmuch as a mathematically rigorous uniqueness 
proof is beyond the scope of this paper, the arguments 
presented will be based either on physical intuition or 
on known criteria of stability of supersonic flow. For 
convenience in referring to Fig. 2 the eight regions 
separated by the lines 17, = 1, A/z = 1, and A2/A, = 1 
will be called regions I to VIII, numbered clockwise 
starting with the triangle OAC. 

(a) Region I(OAC); M, < 1, Mz < 1, Ao/A; < 1.—In 
this region of the chart the ultimate steady flow Mach 
Numbers are subsonic or sonic, the ultimate steady flow 
corresponding to that in a subsonic nozzle for a duct 
with monotonic area change. Inasmuch as an intersec- 
tion of the P;/Q». curve with the A»/A,; curve in this 
region will be accompanied by an intersection in region 
VI, the question arises as to which of the two intersec- 
tions corresponds to the physical solution. If the duct 
has a monotonic area change, the solution in region VI, 
where 1/,; > 1 and A, < 1, indicates that a standing 
shock wave is in the converging section. It is well 
known, from experiments and from the theory of Kan- 
trowitz,? that such a flow configuration is unstable. 
The above argument cannot be used if the channel has 
a centrally located throat, since the shock could possibly 
have a stable position in the diverging portion of the 
channel. The following physical argument, therefore, 
is hypothesized for determination of physical solutions 
for particular initial conditions when the channel area 
varies monotonically or includes a throat. 

This argument is applicable to interactions which 
have initial conditions in regions I and VIII. The 
argument cannot be applied to solutions which occur in 
the remaining regions nor to interactions with initial 
conditions in these remaining regions. The mail 
reason for the former statement in the previous sentence 
is that in the regions other than I and VIII, for most 
interactions, the assumption that the reflected system 
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FINITE AMPLITUDE 


consists Of a W_ wave in the A, section and a |W, 
wave in the A» section (i.e., P; /(Q2 = PQs) is not satis- 
fied. 

The initial steady flow Mach Numbers, .\/; and .\/,, 
which correspond to a point in a certain region, are 
noted on the chart. Assume now that a wave of in- 
fnitesimal amplitude is propagated through the chan- 
nel. If either the initial or final steady flows have no 
shocks in the variable area channel, it appears physi- 
cally reasonable that the ultimate steady flow Mach 
Numbers (.1,, M2) cannot be much different than the 
initial Mach Numbers (\/,, 75); the solution, there- 
fore, must be in the same region as (1/4, 7/5). The pre- 
ceding argument can be extended to waves of finite 
mplitude if these waves may be thought of as a train of 
snall amplitude waves with each wave sufficiently 
separated in distance (or time) so that an equilibrium 
steady flow is attained step by step. 

The use of the principle hypothesized above for the 
letermination of the physical solution for regions I and 
VUI, for cases where both initial and final conditions 
lie in the same region, is as follows: This final conditions 
M, \2) lie on the A2:/A,; = constant curve passing 
through the initial state (73, M.). For example, con- 
sider that the initial Mach Numbers (44, J/g) lie in 
region I and that the incident wave is a downstream 
compression. The point (1/4, A/,) is indicated on Fig. 2 
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by the point a. For an incident downstream compres 
sion wave the ratio P,/Q, will be larger than P;/Q¢s; the 
value of P,/Q» for the assumed wave is denoted by the 
point b on Fig. 2. The point b is therefore the solution 
for (.\,, M2) since it follows the above principle; i.e., 
it lies on the A,/A»2 curve passing through (1/4, 5). In 
creasing the strength of the incident downstream com 
pression further increases P,/Q»2 until for some value of 
strength the point c is reached for which 1/2 is sonic. 
Further increase of the strength of the incident wave 

* will give a value of P;/Q; which lies in region IT; it will 
be shown in paragraph (c) that this value does not neces- 
sarily correspond to P;/Q2. 

Up to now we have assumed that the duct between 
sections 1 and 2 had a monotonic variation of cross sec- 
tion. Consider now that the variation is nonmonotonic, 
and in particular that a throat section exists between 
sections | and 2. Thus if the throat area is the critical 
cross-section area relative to A; (i.e., sonic velocity 
exists in throat section), a solution in region I can repre 
sent Laval nozzle flow with a stationary ‘‘isentropic 
shock”’ downstream of the throat. The stationary isen 
tropic shock can for fixed J/, be anywhere between the 
throat and section 2; it also is unstable in the sense 
that an infinitesimal change of back pressure will dis- 
lodge the ‘‘shock’”’ from the variable area duct. Thus 
todetermine meaningful solutions with stationary shocks 








W4 =Downstream Propogating wave. 
wWi= Upstream Propogating Wave 
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in the variable cross-section duct (initially or finally), 
the real nonisentropic shock relations niust be used, and 
the chart of Fig. 2 is therefore not useful. 

Consider now interaction problems with initial condi- 
tions in either regions II or V with the value of P:/Qs5 
lying in both regions I and VI. The hypothesized rule 
for obtaining the physical solution cannot be used here 
since the initial and final conditions (if the latter are 
given by the chart) lie in different regions. No general 
rule has been formulated for determining the physical 
solution for such interactions. It is generally not dif- 
ficult, however, to determine what appears to be the 
most physically reasonable solution. The appropriate- 
ness of a selected solution may be checked by consider- 
ing the transition to the solution by calculating solu- 
tions for several incident waves ranging in strength from 
zero to the strength of the considered incident wave. 
When a solution is obtained it should be checked for 
violations of boundary conditions (e.g., sonic throat con- 
ditions) or physical conditions such as the laws of wave 
propagation (e.g., a W_ expansion wave cannot move 
upstream in section A; when 1/7; > 1). When such viola- 
tions exist, it is usually a simple matter to select a wave 
scheme which enables satisfaction of all physical and 
boundary conditions. 

(b) Region VIII(OAB); M, < 1, Mz < 1, Az/A; > 
1.—For a duct with monotonic area change, this region 
gives solutions which correspond to flow in a subsonic 
diffuser. It will be noted, however, that for certain 
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values of P,/Q, and A»/ A, there is no intersection of the 
two curves, therefore, no solution given by the chart 
The nature of the solution in these instances may be de. 
termined by the following argument. Consider that the 
initial Mach Numbers (1/4, 1/;) are positioned on the 
chart at point d in Fig. 2. If the incident wave is q 
downstream compression the value of P;/Qs will be 
illus 


ereater than P;/QO;; assume a typical solution 
» - 


trated in Fig. 4a) to be indicated by point e. Increasing 
the strength of the incident compression will increase 
P,/Q» until for a certain strength of the incident wave 
the solution will be at f, at which point 1/7, = 1. 
the A»,/A,; curve has a maximum at this point, further 
increase of P,/Qz» (1.e., the strength of the incident com- 


Since 


pression) yields no intersection of the P;/Q. and A./A, 
curves. A slight increase of the strength of the incident 
compression will now yield a solution (as indicated 
schematically in Fig. 4b) with a standing normal shock 
wave in the diverging channel near station 1, the Mach 
Number at | being unity. Increasing the strength of the 
incident compression will move the standing shock to a 
further downstream position in the diverging channel, 
until at some strength it is stationary at station 2. 
Further increase of the shock strength will cause the 
standing shock to move downstream as an upstream 
propagating wave (see Fig. 4c). Solutions of the above 
type can easily be obtained using the wave schemes of 
Fig. 4. It should be noted that in order to find the 
standing position in the diverging channel, this shock 
will have to be treated as an exact shock, inasmuch as 
The 
solution with real shocks must be obtained by an itera- 


shock position is not fixed for isentropic shocks. 


tive procedure 

When the variable area duct has a nonmonotonic 
change in cross section—in particular a throat—con- 
sideration of the throat boundary condition must be 
given. The discussion of this situation in paragraph 
(a) applies as well here, and will therefore not be re- 
peated. 

(c) Region II(AEDC); M, < 1, Mz > 1, A2/A, <1— 
In this region .V/, is subsonic or sonic, .\/, is supersonic or 
sonic, and A» is less than or equal to A;. If the duct is 
proportioned correctly (i.e., if it has a throat of proper 
area at which the velocity is sonic) the solution here will 
correspond to flow in a Laval nozzle. Otherwise the 
solution indicated corresponds to flow with a stationary 
that the 


If the duct 


expansion shock and we must conclude 
physical solution is not given by the chart. 

has a monotonically decreasing cross section the re- 
quirement for a physical solution is that three reflected 
waves occur (see Fig. 3b); an upstream propagating 
wave W_! moving upstream into the duct of area A), 4 
downstream propagating wave W,* moving dowr- 
stream into the tube of area A», and an upstream propa- 
gating wave W_? which moves downstream into the 
duct of area A». The Mach Number at the throat 
(station 2) is unity, the expansion wave W_? moving 
downstream because of this. The solution can easily be 





obta 
usin 
Me = 
d 
ejthe 
wave 
§ mic 
The | 
f th 
the a 
not § 
nami 
press 
trans 
the ¢ 
shock 
shock 
comp 
until 
streal 
\ sin 
wave 
the 11 
sion 0 
ing sk 
move’ 
inere< 
ward 
into t 
(e) 
alreac 
erally 
gives 
and V 
Fig. 5 


An 
theory 
shock 
area. 

Sine 
ready 
Cussio 
erated 

f hig 
pressic 
very ( 
g10n 0 
measu 
screen: 
peatal 
light Ss 
which 
tern m 
ment. 


1 of the 
chart 
- be de- 
hat the 
on. the 
ve 1S a 
will be 
illus. 
reasing 
icrease 
t wave 
Since 
further 
it com- 
A,/A, 
icident 
licated 
| shock 
» Mach 
1 of the 
*k toa 
lanuel, 
tion 2. 
ise the 
stream 
above 
mes of 
nd the 
shock 
uch as 
The 


1 itera- 


lotonic 
con- 
ust be 
igraph 
be re- 


<1-— 
onic or 
duct is 
proper 
re will 
se the 
onary 
it the 
e duct 
he re- 
flected 
gating 
| Aj, a 
down- 
propa- 
to the 
throat 
10ViNg 
sily be 





FINITE AMPLITUDE 
obtained from this system of waves with ./ determined 
using Eq. (6)) given the area ratio A,/A, and letting 
i= E. 

(d) Region VIII (AIJB); M, > 1, Mz < 1, A2 


1—The steady flow solutions in this region correspond 


A, > 
either to flow in a channel with standing normal shock 
wave or, if the duct has a throat, to shock-free super- 
sonic diffuser flow with sonic velocity at the throat. 
The latter is the only meaningful solution in this region 
of the chart, since it has already been pointed out that 
the assumption of an isentropic stationary shock does 
not give physically useful information. The gas dy- 
namics of the situation is as follows. The incident con 
pression has no upstream propagating reflection, but is 
transmitted and causes formation of a standing shock in 
the diverging portion of the duct. (To calculate the 
shock position, calculations must be made using real 
shock theory.) Increasing the strength of the incident 
compression moves the standing shock further rearward 
until finally for some strength this shock moves down- 
stream as an upstream propagating wave (see Fig. 4d). 
A similar phenomenon would occur with the incident 
wave being an upstream propagating expansion. If 
the incident wave is a downstream propagating expan- 
sion or an upstream propagating compression, the stand- 
ing shock (if we assume one exists initially) would be 
moved forward as the strength of the incident wave is 
increased, until for some strength the shock moves for- 
ward out of the diverging section, propagating upstream 
into the A, section duct. 

(e) Regions III, IV, VI. 
already developed show that these regions do not gen- 
erally give useful flow solutions. The chart, therefore, 
gives physically useful solutions principally in regions I] 
and VIII. 
Fig. 5 for values of y of 1.4 and 1.25. 


V, and The arguments 


An enlarged view of these regions is given in 


EXPERIMENTS 


An experimental test of some aspects of the foregoing 
theory was made by studying the passage of a plane 
shock through a two-dimensional channel of varying 
area. 

Since the theory and operation of the shock tube al- 
teady have been described in detail,'! only a brief dis- 
cussion will be included here. Shock waves are gen- 
trated by breaking a partition which separates regions 
f high and low pressure. When this occurs a com- 
pression wave starts through the low-pressure air and 
very quickly becomes a plane shock followed by a re- 
gion of uniform flow. The velocity of this shock, as 
measured by the interval between pulses from two light 
xreens a known distance apart, is consistently re- 
peatable within one per cent. The pulse from the second 
light screen is also fed into an adjustable delay circuit 
which triggers a short duration spark. The flow pat- 
ttn may thus be observed at any stage in its develop- 
ment. 


WAVE INTERACTIONS 511 





Frevre 5. 


Waraive CHARTS FoR Saaurrvaw a TsauTRoeic 
tle Inter acnon Progiiot Ar Aw Aret DisconTiylFy 


} 





FIGURE 'S. 


A Mach-Zehnder interferometer is used to obtain 
interference fringe patterns from which the density field 
may be evaluated. In two-dimensional flows the fringe 
shift at any point is simply proportional to the density 
change. For most experiments the interferometer is 
adjusted to give straight horizontal fringes. If then a 
picture of the original straight fringe pattern and a flow 
picture are superimposed, those regions where fringes 
cross become readily visible and the entire density field 
can be sketched in quite easily. 
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Fic. 6. 


A special channel for these experiments was made 
from ground stock beveled to provide the change in 
area as shown in Fig. 6. The lower plate could be re- 
versed to provide a diverging channel as well. The 
width of the plates was exactly the width of the tube, 4 
in., so that the sides of the shock tube served as the 
sides of the channel. The maximum depth of the chan- 
nel was 1.01 in. and the ratio of small to larger area was 
0.504. The upstream ends of the plates were beveled 
to a sharp angle to provide a clean entry for the incident 
shock. Experience has shown that under these con- 
ditions the strength of the wave impinging on the area 
change section is then the same as the strength of the 
original shock produced in the tube. 

In order to keep losses from flow separation down to a 
practical level in the diverging case the angle of the area- 
The length Z of 
When seen 


change section was made small, 10°. 
the area-change section was then 2.85 in. 
through the 5-in. windows provided for the inter- 
ferometer this left about 1 in. on the upstream and 
downstream ends of the section for observations of the 
flow into and out of the area-change section. 

In the converging case the corners A and B were 
sharp. In the diverging case corner B at the entrance 
of the area-change section produced wave reflections 
that complicated the flow unnecessarily. For this 
reason it was rounded off to a radius of curvature of 4.5 
in. for the tests with the diverging channel. The time 
available for making observations before the incident 
shock returned from the end of the tube varied with the 
shock strength. In the most severe case, p3/)) = 5.05* 
(i.e., for incident shock pressure ratio) this time was 
about 1,900 microseconds after the incident wave 
reached the area-change section. Examination of the 

* The subscript 0 is used for regions 4 and 6 when the velocity 


is initially zero in these regions. 


Diagram of test channel in shock tube 


interferograms showed that for all shock strengths the 
time available was sufficient for the establishment of 
equilibrium. 

To obtain the proper pressure ratio across the par- 
tition, it was sometimes necessary to partially evacuate 
the portion of the tube containing the test section. The 
initial pressure in the test section is given in Table | 
for each of the tests. The initial air temperature was 
74° + 1°F. for all tests. 

Two series of interferograms were taken at short in- 
tervals showing the transient development of the flow 
from the beginning to the onset of steady conditions. 
This was done with the converging channel at shock 
strengths of 1.96 and 3.95. In addition several inter- 
ferograms were taken after steady flow had been estab- 
lished by shocks of strength p3/po9 = 1.98, 2.52, 3.08, and 
5.05. The times at which the interferograms were 
taken, measured with respect to the time at which the 
shock first reached the area-change section, are included 


in Table 1. 


Selection of Experimental Density for Comparison with 

Theoretical Calculations 

In the theoretical analysis of the wave interactions, 
one-dimensional flow is assumed at any cross section 0! 
the channel. From the density contours in Fig. 15 it 1s 
evident that this condition is not fulfilled in the actual 
case. For a reasonable comparison between equilibrium 
conditions as measured and those obtained theoretically, 
it is necessary to examine this variation of density with 
the purpose of choosing a value of density which most 
closely fulfills the theoretical conditions. The variatiot 
of density measured in these experiments arises mainly 
from the effects of viscosity and heat transfer. The 
density is seen to be greater near the walls than in the 
center of the channel. For a boundary layer with m 
heat transfer, the variation of density would be opposité 
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TABLE | 
Experimental Data 


p 
Pp p Mm. Hg agt /L pPi/p p2/ p | 1 
Converging Channel A./A, = 0.504 
1.96 163 0 1.607 1.00 
0 1.61 1.00 
0.215 1.64 1.00 
0.224 1.66 1.00 
0.473 1.71 1.00 
0.478 1.69 1.00 
0.688 1.68 1.00 
0.716 1.73 ees 
0.950 1.74 Lae 
0.917 1.76 Live 
1.403 1.85 1.75 
1.428 1.85 1.78 
1.897 1.91 1.75 
3.24 1.96 1.75 
4+.70 1.96 L.75 
6.15 2.01 1.79 
7.62 1.97 1.75 
9.24 1.94 1.72 
10.65 1.97 
3.95 65 0.119 2.52 1.00 
0.253 2.58 1.00 
0.358 2.40 1.00 
0.564 2.34 2.80 
0.865 2.89 2.75 
1.296 3.10 2.84 
2.26 3.88 2.78 
3.70 4.07 2.87 
5.16 4.10 2.98 
6.54 4.22 3.04 
7.96 4.32 3.08 
9 41 4.32 3.18 
1.98 161 6.47 1.96 1.73 
8.13 1.96 l.@2 
10.23 1.96 1.74 
oa 91 5.88 2.52 2.08 
7.84 2.49 2.05 
9 85 2.53 2.04 
3.08 154 5.66 3.15 2.46 
7.60 3.17 2.40 
9 60 3.19 2.44 
9.05 32 5.78 5.32 3.69 
7.26 5.46 3.87 
8.73 5.62 3.98 
Diverging Channel 4,/A» = 0.504 
2.49 240 1.90 1.32 1.64 0.54, 0.77 
2.42 1.28 1.64 0.66, 0.72 
2.88 1.2 1.60 0.70 
3.36 | .22 1.55 0.69 
3.87 1.30 1.47 0.69 
3.98 1.35 1.438 0). 682 
4.46 1.33 1.42 0.673 
4.94 1.30 1.43 0.673 
5.36 1.34 1.44 0.668 
6.89 1.29 1.49 0.660 
8.41 1.35 1.51 0.660 
9.47 1.32 1.52 0.651 
2.02 131 1.61 1.28 1.49 
4.18 1.09 1.42 0.525, 0.603 
5.67 1.10 L.35 0.528, 0.588 
7.85 1.14 1.38 0.56 
9.71 1.12 1.38 0.56 
3.40 104 6.21 2 00 1.43 0.89 
7.63 2.04 1.44 0.88 
9 13 2.04 0.87 
1.56 288 6.65 . oF 1.28 
8.67 1.21 1.28 
10.61 1.2] 1.27 
5.05 3107 —0.07 1.00 1.00 
5.66 2.80 0.95 
7.12 2.95 0.90 
8.55 2.90 
to this. With heat transfer, however, the variation of 


density as measured is qualitatively the variation which 


actually exists. This can be shown to be true by calcu- 
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lating the velocity profile from the measured density 
variation. 

The equation used for the calculation of the velocity 
profile was developed by Crocco,'? for the laminar 
boundary layer, and is based on the assumptions of 
zero pressure gradient, linear variation with tempera- 
ture of the viscosity and conductivity coefficients, and 
Prandtl Number equal to 1. The result is also con- 
sidered to hold approximately for a turbulent boundary 
layer. The equation was derived for the case of the 
boundary layer on a plate in a free stream; it is applied 
to the channel in the present case by taking conditions 
at the center of the channel as the free stream condi- 


tions. The equation is as follows: 


—s iC ) 
2 uy 


where the subscripts w and f refer to conditions at the 
wall and in the free stream, respectively. The equation 
may be manipulated to solve for the velocity ratio u /u,: 


4(6,, — 6) 


l — 1 — (9) 
\ (0, + 1)- 


where 
1-— 7/T, | — 7/7; 
6 = , 6, = 
[(y — 1)/2].17,7 [(y — 1)/2]M,? 


With the assumption of constant pressure across the 
section and the perfect gas equation of state, the tem- 
peratures may be related to the density ratios by: 

r Py / po To py/ Po 


T, — p/ po’ f, Pw/ Po 


Thus by means of the measured density ratios and 
with an assumption for the Mach Number at the center 
of the channel, \/,, the velocity profile may be calcu- 
lated using the above equations. 

As an illustration, the method is applied to the data 
obtained with the shock of 3.95 shown in Fig. 13, 
at/L = 9.41. As pointed out above, conditions here 
are very nearly steady. The Mach Number /, used in 
the calculations was that determined from the analysis 
with the exact shock theory. The section for the in- 


vestigation was taken about '/, in. upstream of the 
The variation of density across this sec- 


area change. 
In applying these data to the 


tion is shown in Fig. 7. 
calculation of the velocity profile from Eq. (9), the most 
important quantity is the density ratio of the air at the 
This quantity cannot be determined from 


wall, p../ po. 
A specious method for deter- 


the interferograms alone. 
mining it is to extrapolate the density profile to the wall 
as shown in Fig. 7. It is seen that the velocity profile 
for this assumption resembles that obtained in fully 
developed laminar flow in a channel. Another method 
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of obtaining a value of p,, py is to assume that the tem- 
perature of the air adjacent to the wall is the room tem- 
perature 7). This assumption appears justified, since 
the wall is at room temperature before the test, and can 
increase in temperature but a small amount in the very 
required for the 
process to go to equilibrium. 7) the den- 
sity ratio at the wall was calculated using the theoretical 


short time (less than 2 milliseconds 
Using T, = 


value of pressure ratio p; py found from Table 2 as 
Pw! Po = (Pi/ po) (Tu/To) = 8.1 


[he velocity profile calculated from Eq. (9) using this 
value of density at the wall is also given in Fig. 7. In 
this case the profile obtained could represent either a 
fully developed turbulent flow in the channel or uniform 
flow in the central region of the channel with thin 
boundary layers at the walls. 

Insight into which of the p,,. p, assumptions is correct 
may be obtained by calculation of the possible boundary 
layer thicknesses at the section. The Reynolds Number 
based on length from the channel entrance to the section 
using the data of Tables | and 2) is roughly 2 X 10°. 
This value is below the incompressible transition value 
of 5 X 10°. Since the cooling of the boundary layer 
which takes place would tend to stablize the layer and 
further increase the transition Reynolds Number, the 
existence of a laminar boundary layer appears possible. 
The approximate thickness of laminar boundary layer 
calculated for the test condition is 0.18 in. Since the 
channel height is 1.01 in., the assumption of extra- 
polated p/p», which indicated fully developed laminar 
channel flow, is definitely ruled out. The possible ex- 
istence of a turbulent boundary layer cannot com- 
pletely be denied, since disturbances which might pro- 
duce early transition may be present. (The height of a 
completely turbulent boundary layer was calculated as 
0.51 in.) In any event the assumption of room tem- 
perature for the air adjacent to the wall is consistent 
with velocity profiles which appear possible, thus the 
variation of density across the channel at “equilibrium” 
times is attributed mainly to heat transfer rather than 
to variation of velocity. Clearly the fluid at the center 
of the channel will be least affected by heat flow out 
through the channel walls and this is closest to the 
theoretical condition of a completely insulated channel. 
For this reason values of density measured at the center 
of the channel are used for comparison with the theoreti- 


cal results. 


Corrections to Density Measurements 

The determination of density from an interferogram 
is predicated on uniform conditions along the path fol- 
lowed by each light ray crossing the channel. Actually 
any ray must pass through the boundary layers on the 
observation windows and may also be bent by density 
gradients perpendicular to the light path. Since 
these effects are small they have been treated separately. 
Calculations of the refraction error using the theory of 
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Wachtell'® and Howes and Buchele' show that refrac 
tion may be safely ignored in all regions except the im 
mediate vicinity of the channel walls. The boundary 
layers on the windows increase the apparent density 
about 3 per cent in the latest pictures. No correction 


was made for this effect either. 


Time Variation of Density 

Examination of the transient flow conditions pro 
duced by wave interactions at the change of area sec- 
tion allows the determination of the times required for 
the establishment of steady flow conditions. The ap- 
proach and attainment of steadiness in the converging 
channel are illustrated by the plots of density ratio 
versus time (measured in units of dof ZL) after the im- 
pingement of the initial shock wave on the area-change 
section, shown in Figs. Sand 9. These density ratios are 
averages taken across the channel (except when the 
heat transfer effect is noticeable) a quarter inch up 
stream and downstream of points A and B in Fig. 6, 
respectively. It is seen that steadiness is practically 


obtained at an dof L value of about 3 
The rise in density 


approximately 
650 microseconds) for both cases. 
after the apparent attainment of steadiness in Fig. 9 is 
attributed to the heat transfer effect discussed pre- 
vil nusly. 

A theoretical time-density variation is also shown in 
Fig. S (initial shock pressure ratio is 1.96) for com- 
parison with the experimental results. To arrive at the 
theoretical results of this section, a method of character- 
istics solution was used to solve the one-dimensional 
unsteady equations of motion in a channel of the same 
shape as the experimental model. This method is 
similar to that discussed by Guderley' in which all 
waves are considered to be isentropic; however, a re- 
finement to Guderley’s treatment was employed allow- 
ing the initial wave to be treated as an adiabatic normal 
shock as it moves through the channel while all other 
waves and flow regions are still considered to be isen- 
tropic.’ Investigation of the actual waves in the con- 
verging section from the test interferograms showed this 

The resulting 
shock 


to be a desirable feature of the analysis. 


network for an initial pressure 


characteristic 
ratio of 1.96 is shown in Fig. 10. 

Density ratios were calculated varying with dot/L at 
the upstream and downstream ends of the converging 
section by using the isentropic relations behind the 
transmitted shock waves and the normal shock relations 
across it. It is seen (Fig. 8) that at the upstream end 
of the channel the theoretical and experimental agree- 
ment is good except at the upper limit of the transition 
flow region. The experimental values for the down- 
stream or smaller end of the channel also show good 
agreement with the theory. 

Fig. 11 compares experimental and theoretical varia- 
tions of density along the channel at various times dur- 
ing the transition period. The experimental points are 
averages taken across the channel at stations chosen at 
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equal area ratio increments. It is seen that there is re- 
markably good agreement in the shape of the curves. 
The formation of steady-state values upstream through 
the channel is clearly a function of time, the actual ex- 
perimental times being almost identical with those pre- 
dicted by the one-dimensional calculation. 

The observed experimental and calculated positions 
of the initial shock in the converging channel as func- 
tions of time are shown in Fig. 12. The experimental 
points lie consistently to the right of the theoretical 
curve. Since the error is fairly constant for most points 
and since the scale is greatly magnified (the maximum 
error corresponds to approximately 15 microseconds), it 
is believed that the discrepancies are caused by errors in 
the zero time determination and other experimental in- 
accuracies. 

From the one-dimensional calculation, it is interesting 
to note that the W+ waves, after having been trans- 
mitted downstream through the channel, remain essen- 
tially parallel behind the transmitted shock in the time- 
distance plane (Fig. 10). This indicates that the 
transmitted shock is neither strengthened nor weakened 
appreciably after passing through the converging por- 
tion of the channel. Experimentally observed points 
also have this property of essential steadiness at the 
downstream end of the channel immediately after the 
transmission of the initial shock wave; not only for the 
initial shock pressure ratio corresponding to the calcu- 
lation value (Fig. 8), but for the higher shock strength 
as well (Fig. 9). A comparison of Figs. 8 and 11 also 
indicates that practical steadiness is achieved at the up- 
stream end of the channel shortly after the last IW. 
wave reflected upstream from the initial shock, while 
it is still in the converging section, reaches the up- 


stream end. 


Two-Dimensional Wave Interaction Phenomena— 

Discussion of Interferograms 

The two-dimensional reflection of shock waves from 
an inclined surface is an interesting phenomenon which 
has been studied extensively in the Princeton Shock 
Wave Laboratory and elsewhere. A comprehensive 
treatment of the theory and many experimental data 
have been given by Fletcher, Taub, and Bleakney.'* 
In this case the reflections are confined between the top 
and bottom of the channel so that one can follow the 
repeated reflections of the same wave, the formation of 
the composite reflected wave, and the strengthening of 
the transmitted wave. Since these investigations are 
not the main concern of this paper, however, they will be 
treated only briefly. 

Fig. 13 is a series of density contour drawings made 
from test interferograms for flow in the converging chan- 
nel at various times and at the same initial shock pres- 
sure strength. The first drawing shows the shock just 
entering the converging channel and the beginning of a 
Mach reflection at the sloping wall. It will be noticed 
that the cylindrical reflected wave and the density con- 
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a. tours inside it in picture 2 are almost an exact magnifica- bottom of the channel. The formation of the upstream 


a tion with time of picture | in agreement with the theory — reflected wave can be traced through the series of pic- 
> ff Mach reflection.'® tures quite clearly; picture 3 shows the steepening of 
The third picture shows the reflection of the cylindri- the initial cylindrical reflected wave as it moves up- 
cal wave from the top of the channel, and in picture 4, stream; picture 4 shows the increased steepening of the 


this wave has just started to reflect again from the initial reflected wave and the onset of upstream portion 
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of the first wave reflected from the top of the channel, 
as well as the near normalcy of the transmitted wave as 
it leaves the channel; picture 5 shows the gradual con- 
vergence of these reflected waves, the onset of the up- 
stream portion of the wave reflected from the bottom of 
the channel (seen in picture 4), and the formation of 
another upstream moving wave; picture 6 shows further 
reflection in the converging section and the tendency of 
these waves toward a normal attitude as they move up- 
stream; picture 7 shows the reflections in the converg- 
ing section dying out with the last few waves moving 
upstream at the large end of the channel; picture 8 


shows essential steady flow established in the entire 


channel; pictures 9, 10, 11 show the increasing effect 
with time of the heat transfer through the top and bot- 
tom of the channel; picture 12 shows the initial wave 
returning into the test section after having been re- 
flected from the downstream end of the model. The 
actual coalescing of the upstream moving compression 
waves is not actually visible in any of these pictures. 
However, their tendency to become normal is obvious, 
and it is known that a compression wave moving in a 
straight channel will always overtake 

Therefore, 


another wave 
moving in the same direction. the formation 
of a final normal reflected shock wave upstream of the 
area-change section seems to be absolutely predicted by 
these experiments. 
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TABLE 2 
Theoretical Results for the Converging Channel 


§22 JOURNAL OF THE 
Isentropic Wave Theory 

Wave Scheme— a a a a 
p/P 1.98 2.52 3.08 3.95 
M 0.47 0.62 0.74 0.89 
pi/p 1.33 1.55 1.74 2.05 
ps/P 2.24 2.96 3.60 $.55 
M, 0.24 0.29 0.51 0.32 
p1/ po 2.00 2.65 3.30 4.45 
Me 0.55 0.72 0.84 0.98 
p2/ Po 2.24 2.96 3.60 4.55 
p2/ po 1.78 2.18 2.50 2.94 


Note: Wave scheme refers to diagram in Fig. 3 with a region of no flow, represented by the symbol 0, replacing ri 


TABLE 3 
Theoretical Results for the Diverging Channel 


Isentropic Wave Theory 


Wave Scheme—> a bh Cc 

p3/ po 1.56 2.03 3.40 
M; 0.3 0.48 0.80 
pi/pP 0.81 0.63 0.79 
Ps/ Po 2.91 
be/ Pi 

M, 0.48 1.00 1.00 
pi/ pe i Y 1.06 2.02 
M, 0.22 0.31 2.19 
p2/ Po 1.37 1.58 0.48 
p»/ po 1.25 1.38 0.59 


Note: Wave scheme refers to diagram in Fig. 4 with a region of no flow, represented by the symbol 0, replacing regi 


It is interesting to note that there are practically no 
observable reflections taking place in the downstream 
portion of the converging section throughout the tran- 
sient period in the channel. In fact, immediately upon 
passage of the transmitted shock, the flow in this sec- 
tion becomes almost one dimensional and steady while 
unsteady wave reflections continued vigorously through 


the rest of the channel. 


Steady-State Conditions 


Calculations were made to determine the final steady- 
state conditions existing in both a converging and a 
diverging channel (each with a ratio of smaller to 
larger area of 0.504) which are initially under no flow 
conditions and which are impinged upon by normal 
compression waves. The strength of the initial wave in 
both cases varies from | to 6, therefore including all the 
wave schemes of Figs. 3 and 4. The calculations were 
made for both isentropic compression waves and shock 
waves. For the case of isentropic waves in scheme a of 
both Figs. 3 and 4, the method of calculation used was 
that developed in the theoretical section utilizing Fig. 
5a. For all other cases, trial and error methods were 
used.!® Calculated results for initial shock pressure 
ratios corresponding to test conditions are given in 
Tables 2 and 3. Since density is the most easily ob- 
tained experimental property, calculated steady-state 
density ratios are plotted in Figs. 14 and 15 as fune- 
tions of initial shock pressure ratio for comparison with 
experimental results. 

The experimental points given in Figs. 14 and 15 are 
taken from the interferograms. Each point is an aver- 


Shock Wave Theory 

h a a a a h 
5.05 1.98 2.53 3.08 3.95 5.05 
1.03 0.47 0.62 0.74 0.89 1.02 
2.44 1.23 1.54 1.74 2.05 2.40 
6.05 2.26 ».94 3.64 bon 5.95 
0.22 0.24 0.28 0.30 0.31 0.32 
6.05 1.99 2.57 3.19 $ 11 5.19 
1.00 0.55 0.71 0.84 0.97 1.00 
7.55 2.26 2.94 3.64 1 7 6.84 
4.07 Lae 2.33 2.40 2.79 3.45 

gions + and 6 
Shock Wave Theor 
d a b c d 

5.04 1.56 2 03 > +0) 5.04 

1.03 0.3 0.48 O80 1.03 

1.00 0.76 0.53 0.78 1.00 

+. 90 2.ae 3.81 

> SI) 

1.03 0.52 1.00 1.00 1.03 

So. 4a 1.13 1.04 1.92 2.83 

2 19 0). 23 0.35 2.19 2.19 

0.94 1.38 1.69 0.48 0.94 

0.96 1.25 1.39 0.57 0.85 


t and 6. 


age value of the density ratio across the section indi- 
cated in the note accompanying the figure until the heat 
transfer effect becomes noticeable, whereupon the value 
at the center of the channel was used. The number as- 
sociated with each point is the time at which the inter- 
» to the time 
of impingement of the initial wave upon the 
In these figures, only interferograms in 


ferogram for that point was taken relati' 
area- 
change section. 
which steady state appears to have been established are 
used to obtain average point values. At most initial 
shock strengths, only a few pictures at later times were 
taken, the criteria for steadiness being the absence of 
unsteady waves and the reproduction of similar density 
patterns, within heat transfer and boundary layer 
growth limits, over a sequence of three pictures taken 
The upstream and down- 


ible 1 for 


at reasonable time intervals. 
stream average density values are listed i: 
all experimental runs. 

Referring to Fig. 14 it is seen. that, for the converging 
channel, experimental results agree wel! with the 
theory, particularly at low initial shock strengths. In 
the higher range of shock strengths, the heat transfer 
effect is seen; the density value increases 
slightly but continuously with time, even at the center 
The discontinuity in the 


“steady” 
of the channel. slope of the 
theoretical curves at a p;/po value of 4.05 
the change from wave scheme @ to wave 


is caused by 
scheme D in 
Fig. 3. 

Experimental densities in the upstream portion of the 
area-change section of the diverging channel (Fig. 15a 
also agree well with the theoretical predictions. The 
points have slightly but consistently higher values than 
those calculated by the one-dimensional! shock wave 
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theory throughout wave schemes 4, b, and ¢, thus indi- 
cating that the actual reflected expansion wave (only 
present in these three wave schemes) is slightly weaker 
than that predicted theoretically. The downward con- 
cavity of the curves in wave scheme @ is easily ex- 
plained; the initial positive slope is caused by the 
strengthening of the initial compression wave, and the 
continuously decreasing slope thereafter is caused by the 
increasing strength of the reflected expansion wave 
which eventually causes the density ratio, p;/ po, to drop 
below | as \/; approaches 1. The absence of the heat 
transfer effect (there is no consistent increase of p/ po 
with time) is due to the fact that the reflected wave in 
the diverging channel is always an expansion wave 
which keeps the temperature level low relative to runs 
at the same initial shock pressure strength in the con- 
lables 2 and 3). 


verging channel (see “ 

Comparison of experimental and theoretical resuits 
for the downstream portion of the diverging channel 
(Fig. 15b) again shows good agreement except at initial 
shock strengths near the transition between schemes b 
and c. The deviation in this range is undoubtedly 
caused by the existence of boundary layers on the chan- 
nel walls which effectively change the area ratio of the 
channel and by the departure from one-dimensionality 
of the stationary shock wave existing in the diverging 
section as the initial shock strength, );//o, increases. 
Fig. 16 shows interferograms of the steady shock pattern 
for two different incident shocks. 


An interesting plot, Fig. 17, shows the difference be- 
tween predicted (one-dimensional shock wave theory) 
and observed locations of the stationary shock located 
in the diverging channel as a function of initial shock 
pressure ratio. The experimental location is deter- 
mined by the position of the central, or normal, position 
of the shock pattern in the area-change section. It is 
believed that the difference between experimental and 
theoretical results is caused again by the variance in the 
actual flow from the one-dimensional theoretical model. 
The boundary layer build-up with time in the down- 
stream section of the channel, which will move the 
shock wave upstream slightly, is noticeable in the times 


of the various test points. 


CONCLUSIONS 


A theoretical and experimental investigation has been 
made of wave interactions at a change of area in a duct. 
A new method is described by which the steady-state 
conditions after passage of an isentropic wave may be 
determined from a chart. Possible isentropic flows are 
discussed along with arguments for the existence of 
only certain solutions based on physical considerations. 
The region for which 4/7, and M2 < 1 is shown to contain 
most of the physically useful solutions. Adaption of 
the method to other flow conditions is made by in- 
cluding the properties of shock waves in the calculations 
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and by allowing for the possibility of standing shocks in 
a diverging section. 

A characteristics-type solution was used to solve the 
one-dimensional transient flow equations for the case of 
a shock of pressure ratio of 1.96 incident on a converging 
section of area ratio 0.504. The incident wave was 
treated as an adiabatic normal shock. 
flow was assumed to be isentropic. 

Experiments were carried out in the shock tube to 
test the foregoing theoretical developments. A two- 
dimensional channel of area ratio 0.504 was made in 
such a way that the wave could be incident on either 
Density fields were measured with the aid of a 


The rest of the 


end. 
Mach-Zehnder interferometer. 

A series of interferograms taken during various 
stages of the flow show the transient solution to be in 
excellent agreement with experiment. Two-dimensional 
wave effects were small since the angle of convergence 
was only 10°. Viscosity and heat transfer to the chan- 
nel walls formed an appreciable boundary layer in the 
later pictures but theoretical considerations and _ the 
observed steady-state flow fields indicate that densities 
observed in the center of the channel could be properly 
compared with the one-dimensional theory. 

The establishment of steady flow was observed for a 
wide variety of shock strengths using both the converg- 
ing and diverging models. The isentropic wave theory 
fits the measured densities well for weak shocks. For 
strong shocks satisfactory agreement is also obtained 
when the assumption of constant entropy is dropped. 
The usefulness of the new theory for predicting the 
steady flow conditions after the interaction of a wave 
with a change-of-area has thus been verified for shocks 
of moderate strength incident on sections of area ratio 


up to 2:1. 
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Nonsinusoidal Buckling Modes of Sandwich 
Plates’ 


J. N. GOODIER? anp C. S. HSU? 
Stanford University 


ABSTRACT 


Both wrinkling and ‘‘quasi-Euler’’ buckling modes of sandwich 
plates have hitherto been taken as sinusoidal. In this report the 
possibility of nonsinusoidal modes is investigated 
material considered is the special orthotropic one which does not 
resist deformation in its own plane, but only extension and shear 
This choice is made for two 


The core slab 


ing perpendicular to its plane. 
reasons: (1) such an idealized core material has been adopted by 
previous investigators as closer to certain real cores than the iso- 
tropic material, and also because the ‘‘in the plane’’ stress com- 
ponents are believed to play but a minor part in the core action, 
2) it permits a relatively simple analysis of nonsinusoidal defor- 
mations and the satisfaction of end conditions which would pre- 
sent formidable analytical difficulties if the core were isotropic 
It is found that the critical wrinkling loads corresponding to 
sinusoidal modes are not the lowest consistent with ‘‘pinned”’ 
Loads about half as great correspond to buckling modes in 
The main 


ends 
which the deformation is confined to the end zones. 
body of the panel remains flat but undergoes a lateral displace- 
ment. The effect of clamping the ends is also examined 

The pattern of the theory follows closely, though not exactly, 
that of the long column with elastic foundation considering end 
conditions other than simple support. Ratzersdorfer‘ showed 
that when the ends are allowed to deflect the critical load is 
halved, and the deflection localized at the ends. The problem is 
the same when the ends are pinned but the rigid base of the elastic 
foundation is freed to ‘‘float’’ under no forces other than those of 
the springs representing the foundation. This is a close analog of 
the sandwich problem. 


(1) INTRODUCTION 
, ee ASSUMPTION OF sinusoidal modes of buckling in 
problems of elastic stability with respect to in- 
fnitesimal deformations is almost as commonplace as 
the assumption of sinusoidal (in time) normal modes in 
vibrating systems with linear elasticities. Whereas the 
latter is justified by a theorem in the general theory of 
such vibrations, the former remains an assumption 
unless we obtain a general solution of the equations of 
the particular problem and thus, or otherwise, are able 
to show that only sinusoidal modes can occur. 


Received July 28, 1953. Revised and received February 18, 
154 

* The investigation reported in this paper was sponsored by the 
fice of Naval Research, under Contract N6-ONR-251 with Stan- 
ford University (Project NR-064-241) 

1 Professor of Engineering Mechanics 

t Research Assistant in Engineering Mechanics. 
Engineer, International Business Machines Corporation, Pough- 


Keepsie, N.Y. 


Now Research 


A notable example of this assumption occurs in the 
stability problem of the sandwich plate. Ordinary 
theory of elasticity is used to express the behavior of the 
core slab, ordinary thin plate equations for the thin 
metal face sheets. The resulting set of equations is 
intractable without the assumption of, or equivalently 
the limitation to, sinusoidal modes. Naturally this 
leaves open the possibility of nonsinusoidal modes at 
lower critical loads. In this paper this possibility is 
explored by means of a simplified core slab theory which 
has been introduced for other reasons by Reissner,! 
Hemp,’ and Hoff and Mautner.* Nonsinusoidal wrin- 
kling modes are found for simple compression of the order 
of half the critical loads for sinusoidal modes. These 
results may account for some apparently unduly low 
buckling loads which have been found in tests. They 
show that the sinusoidal assumption is not always justi- 
fied in the sandwich wrinkling problem, and that the 
relation of formulas based on it to tests requires re- 
examination. 

Another notable example of the sinusoidal assump- 
tion occurs in many accounts of one of the basic ele- 
mentary problems—the compressed bar on an elastic 
foundation. When the overall length is infinite there is 
a finite wave length of sinusoidal buckling at a finite 
critical load. But, 6 
shows, there are also nonsinusoidal modes (correspond- 
ing to diferent end conditions at infinity) with critical 
loads as low as one-half that of the sinusoidal mode. 
and modes 


as exhaustive investigation* ® 


Fig. 1 indicates the critical loads of 
buckling for (a) pinned ends, (b) clamped ends, and (c 
free ends. In the last the ends are free to deflect, the 
end forces remaining horizontal. The base of the foun- 
dation is fixed. If it is free except for attachment to 
the foundation material, and also rigid and weightless, 
the ends of the column may be pinned without chang 
ing the problem (Fig. 1d). Fig. le shows another foun- 
dation and column underneath the same baseboard, the 
deflections of the column being reversed with respect to 
this foundation. The combination forms a simplified 
model of a sandwich with simply supported face sheet 
edges. The baseboard in the middle represents the 
middle undeformed zone of the core slab in a wrinkling 
displacement. But here the wrinkling is localized at the 
ends, and the baseboard, or the idle middle zone, has a 
vertical rigid body displacement with respect to the 


pinned ends. 
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Fic. 1. Modes of buckling for a column on an elastic founda- 
tion: (a) pinned ends; (b) clamped ends; (c) free ends; (d) ends 
pinned, rigid foundation base floating; and (e) case (d) doubled 
as a sandwich model. 










































(2) THe SANDWICH PLATE—SIMPLIFIED CORE THEORY 


When the core slab of the sandwich plate is isotropic 
and is treated by the ordinary equations of isotropic 
and linear elasticity, the analysis is cumbersome. When 
the buckling deformation is not sinsuoidal the difficul- 
On the other hand actual 
The honeycomb 


ties are all but prohibitive. 
core materials are often not isotropic. 
type, and balsa wood with grain normal to the faces of 
the slab, are easily deformed in the plane of the slab,* 
the constraint exerted on the metal face sheets being 
probably attributable in the main to the elasticities 
expressed by the Young’s modulus normal to the faces 
(Z.) and the shear modulus for shearing out of the 
original plane (G,). References 1, 2, and 3 deal with 
bending and buckling of the sandwich plate on this 
supposition. We adopt it here, and show that with 
such a core it is possible to have nonsinusoidal modes of 
buckling at loads substantially lower than those of the 
sinusoidal modes. The simplicity of the theory for this 
kind of core follows from the possibility of immediate 
general integration of the differential equations of the 
core slab, and of putting the results in the form of 
simple constraint equations for the face sheets. 

The stiffness of the core with respect to loading in its 
own plane being negligible we have, with x-y-axes in 


* It has been pointed out by W. Fliigge that a regular honey- 
comb slab is stiff with respect to equal principal stresses in the 
plane of the slab. We have not found it possible to take this stiff- 
ness into account without losing the essential simplicity of the 


theory given here. 
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the mid-plane of the slab, and the z-axis perpendicular 
to it 


Or = Oy = Try = O (] 


og, and o, being normal components of stress, 7.,, the 
shear stress component in the plane. The conditions 
(1) define what L. N. G. Filonf called ‘‘antiplane stress’ 
(plane stress being defined by a, = r,, = 72, = OQ), and 
it will be convenient to refer to the present type of core 
as an “‘antiplane core.’’ The remaining three stress 
components are T;,, Tyz, ¢, and we write 7., 7, now in- 


stead of 7,2, r,.. They satisfy the equilibrium condi- 


tions 

Or, 0 Or, 0 Or, 4 Or, 4 Oo: 0 , 
= ; = ‘ = (7) 

Oz Oz Ox Oy Oz 


With u, v, w as the displacement components, the stress- 
strain relations are 


—— (= 4 ~*) 
Ox dz /’ 

_ , (ow *) ,_ wi. 

7, = G, (= + dz)” og; = FE, - (3) 


Eqs. (2) show that o, must be linear in z, so if a; and a; 
be the values of o, at the upper and lower faces of the 


slab, z = +h/2, we have 


oz = (1/2)(0, + o2) + (l/h) (oO, — on)z (4) 


With this the last of Eqs. (2) becomes 


Or, OT, l : 

mas + ie + (o, — o2.) = O (5) 

Since 7., 7, are, by Eqs. (2), independent of z, they 
do not vary through the thickness, and may be regarded 
as the values of the shearing stress components applied 
to the faces of the slab. These applied shears must then 
be equal at corresponding points of the two faces. 
There is of course no such requirement for a slab of 
ordinary isotropic material, its loading being subject 
only to the requirement of overall equilibrium. Now 
Eq. (5) is a relation between the loadings 7., 7,, 01, 72 0 
the two faces, and therefore the ‘‘antiplane’’ core slab 
can withstand only such face loads as satisfy it. Again 
there is no corresponding condition for a slab of iso- 
The only admissible edge loading for 
vertical shears on 
with 


tropic material. 
the antiplane slab consists of r., 7, 
the edges of the horizontal slab 

mined by the edge values of 7., 7, as shears applied to 


values deter- 


the faces. 
We next integrate Eqs. (3) for the displacements. 


The last one, using Eq. (4), gives 
2Ew = (0, + oo) + (1/h) (0; — oo)2? + 2E.w (6 


where w is an arbitrary function of x and y. Using Eq. 


(6) in the first two of Eqs. (3), recalling that 7., 7, are 


+ L. N. G. Filon, Proc. Roy. Soc., London, Ser. A., Vol 


137, 1937. 
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independent of z by Eqs. (2), we find 


Tz TW ] O 
Eu = 2k ( _ s— —s* — (oe: + oe) — 
= 4 <f a 


. Ay Ty OW l . O 
9k yy = ZE : — ) = = (o; + G2) — 
G, oy 2 Oy 


0,1 — ov) + 2Ew, (8) 


3h Oy : 
where % and vw are arbitrary functions of x and y. 
With Eqs. (6), (7), (8), and (5), 7, and 7, being inde- 
pendent of 2, all the differential equations, and the 
boundary conditions, are satisfied. The displacements 
lue to face loading 7., 7,, 01, o2, and the edge shear 
loading implied by r., r, are determined by Eqs. (6), 
7), and (8) except for the arbitrary functions (of x and 
y) Mo, Y, and wo. This arbitrariness follows from the 
ibsence of resistance to deformations in the plane of the 
core—the core can for instance be stretched in its plane 
by zero forces. The general deformation which the core 
admits without load is obtained by putting 0, 0, 7., Ty 
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eliminating the arbitrary parts expressed by wu, v%, and 


Wo. 


Thus Eqs. (6), (7), and (8) yield 


2B. ut = £(1/2)h (or + 02) + 
(1/4)h(o, — oo) + 2E aw (6a 
. f 41 ‘ Tz OW l 4 Oo 
2E = +E h nT h? (0, T 02) F 
ts G Ox S Ox 
|... @ 
h* (ao; — oo) + ZEA (7a) 
24 Ox 
jv Ty OW l m re) 
»]} = I h . _ —_ h? is + @ + 
V2 G, Oy S Vv 
1 0 
h? (0, — oo) + 2Eaw (Sa) 
24. Oy 


From these six equations the arbitrary functions can be 
Formation of 4, — us and 2; 
The derivatives of wo are found from Eq. 


eliminated. — v removes 
Uy and %. 
(6a) by adding these two equations and differentiating. 
Finally wo is eliminated from Eq. (6a) by subtraction. 
We find then the following three equations relating the 
face displacements of the antiplane core slab to the loads 


causing them 











equal to zero in Eqs. (6), (7), and (8). Then a h 
2(u; — U2) +h (wm + Ww) = 2-7. + 
w = Wh, u = —32(Ow)/Ox) + uM, Ov G 
y = —2(dw/dy) + tw (9) Pass 
° 3(Owwe/Oy) + % —— ae (a, — ao») (10) 
; 6E 
These have the same form as the displacements of a 
plate in ordinary thin plate theory, lines normal to the h 
J . e 7 ; i 3 ‘ 2(v, — vo) th (w, + Wo) = 2-7, + 
middle surface remaining so in the deformation. ! G 
In spite of this free deformability the antiplane core is | ra) 
able to exert a stiffening restraint on the face sheets 61 h’ he (0, — a2) (II) 
bonded to it, and this restraint can be expressed by ‘ 
evaluating the displacements for z = +//2, and then 2E(w, — We) = hlo, + o (12 
4 N 
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Rewriting Eq. (5) as 
Or: OT, l 
(ao, — oo) = O (13) 


O. oy h ij 


we have here a set of four equations of constraint to 
apply to the face sheets. The isotropic core slab would 
not of course yield any such point relations between its 
boundary forces and displacements. 

When the deformation is symmetrical with respect to 
the mid-plane z = 0, we shall have m4 = um, v1, = 2, 
If the normal loading on the faces is also 


W, = —We. 
symmetrical, o, = 2, and it follows from Eqs. (10) and 
(11) that 7, = r, = 0. Eq. (12) becomes the only sur- 


viving equation of constraint, and it reduces to 


a, = (2E./h)w, (14) 
Thus in this case the antiplane core supports the face 
sheets exactly as the ‘elastic foundation” of beam and 
plate* theory, giving a supporting reaction per unit 
length or area proportional to the deflection, and the re- 
sults of this theory can be taken over. 

For instance each of the modes of buckling and criti- 
cal loads of the bar on the elastic foundation has a 
counterpart in the symmetrical buckling of the rectangu- 
lar sandwich plate under simple compression. The bar 
corresponds to a strip of face sheet of unit width. P is 
the load per unit width in one sheet, and the foundation 
modulus isk = 2/./h. Instead of £/ we shall have D = 
E ,t®/12(1 — v/*) where tis the face sheet thickness, /; its 
Young’s modulus, and v, its Poisson’s ratio, It is as- 
sumed that no anticlastic curvature occurs. 

Thus corresponding to the modes and critical loads of 
the bar indicated in Figs. la, b, c we have the sym- 
metrical modes and critical loads (per unit width for 
one face sheet) indicated in Figs. 2a, b,c. The values of 
the latter are expressed directly as P.,, and also in 
terms of the dimensionless quantities 


Ps Bs £E : 
<= oS *- 
th, G, h 


G. h m E 
~ . where Lk, = f Z 
Ek, t ] ae 


(15) 
which are employed again later. The question remains 
as to which, if any, of these modes is likely to occur in a 
test. Fig. 2c can be ruled out since it requires a spread- 
ing apart of the ends of the face sheets, and this in all 
probability would be prevented either by a framing 
around the sandwich panel or by frictional contact with 
the compressing heads of the testing machine. Only an 
infinitesimal amount of frictional force would be ade- 
quate to convert Fig. 2c to Fig. 2a or b, for an infinitesi- 
mal amount of buckling. Thus Figs. 2a and b represent 
the probable modes, and they have the same critical 
load (for large length to wave-length ratio as in wrin- 


kling modes). We shall, however, find lower loads asso- 


* See for instance S. Timoshenko, Theory of Plates and Shells, 


pp. 30 and 248, 1940 
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ciated with the antisymmetrical modes which we now 


consider. 


Mopes ANTISYMMETRICAL ABOUT THI 
MIDDLE PLANE 


(3) SANDWICH 


Again the plate is uniformly compressed in the length 
wise (x) direction, and anticlastic curvature is supposed 
absent. Buckling then consists of a plane deformation, 
the displacement v being zero, and u and w being inde- 


pendent of y. The stress components in core slab and 


face sheets are also independent of y, and 7, = 0 by 
Eqs. (3). 
The core Eqs. (10), (11), (12), and (13) become 
| d 2 
2(u, — we) +h (Ww, + Wo) = hr, + 
dx G, 
h? d 
ae (0, — G2 16) 
6, dx 
2E (W, — We) = h(o; + 02 17 
dr,/dx = — (1/h)(o, — o2) (18) 
But the restriction of the present considerations to 
modes antisymmetrical about the mid-plane z = 0 
means 
“le = —), We = Wi, 
and thus from Eq. (17) 0. = —o,. So Eqs. (16) and 
(1S) reduce to 
2u, + h(dw,/dx) = (h/G.)r. + (h?/6E.)(do,/ dx) (19) 
dr,/dx = —(2/h)o (20 


these now being the only equations of constraint to be 
applied to the face sheets. 

The equations of the face sheets have the form of 
beam-column equations and are (Fig. 3 shows the con- 
dition of an element of the upper sheet) 


(dM/dx) — Q — 7.(t/2) = 0 (21) 
dQ | dw 
+o, + (P — pit) —- = 22 
dx dx? 
M = D(d*w/dx?) (23 


Here ./ and Q are bending moment and shearing force 


per unit width. The term p,f in Eq. (22) expresses the 


m,, @ 
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change of thrust P which appears on buckling, and is 
We may however dis- 


required for equilibrium with r,. 
Then we have by 


eard it for infinitesimal buckling. 
elimination of Q and \V/ 

d*w t dr, 
> 


TO = 0 


d o/h) 
) — z 
dx4 dx? 2 dx 


I (24) 

For the sake of simplicity it is preferable to restrict 
the present analysis to very thin face sheets, for which 
we can take* u, = 0, neglect 7.(t/2) in Eqs. (21) and 
24), and also, of course, identify the deflection w of the 
sheet with the displacement w, of the core. We have 
then from Eqs. (19), (20), and (24) 


(dw/dx) — (7./G.) — (h/6E.)(do,;dx) = 0 (25) 
(dr,/dx) + (2/h)o, = 0 (26) 
D(d*w/dx*) + P(d?w/dx?) + o, = 0 (27) 


These are three equations for w, o, rz. Eliminating o 


and r, we have for w itself the equation 


d*w ( P 
dx® D 


IZE ) d toy 
h?G, dx* 
12K, ( hG, P ) d*w 


: — (28) 
h°G, dx? 


= 0 
2D D 
Once w is known, o; and +, follow from Eqs. (27) and 
26). The solutions of Eq. (28) as an equation in 


19 9 ° ° Ax e 
d’w/dx* are of the form e“™ where \ may have the four 
values determined by 


iis os 12\2?/F | 
+ Bisons (2) (£- 5s) a0 
t2 t? 2 


the coefficients of Eq. (28) having been replaced by the 
dimensionless quantities, Eqs. (15). 
The four values are 


Mu =1V6/t VS — C)- V(S+ C)? — 2CF =r 


(29) 


V (S + C)* — 2CF = p 
(30) 


1vV6/t YS — C) + 
Ny = =-_ 
Then the general solution of Eq. (28) may be written 


cosh rx + A3 sinh px + 
A, cosh px + Asx + Ag 


w= A, sinh rx + A» 
(31 ) 
the constants A;, A», A;, and A, being complex as re- 
quired to make w real. 

The character of the deflection w evidently depends 
on whether r and p are real, imaginary, or complex. 
We now construct a simple diagram which shows these 
features at a glance when S, C, and F are given. 

From Eqs. (29) and (30) we can write 
t yr? : ; 

' =1l1—s+V(1+s)’?-4f (32) 
6C (2 
* This is accepted simplification for wrinkling modes, but not 


lor overall, quasi-Euler, modes 
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Fic. 4. Character of solutions of Eq. (28 
where 
S P.,G hr F Gh? 
‘= = ——_., f= : ——- 
. E,E ag : 2¢ 2E Et 


It is useful to consider points in an f, s-plane, both f 
and s being positive, and to draw the parabola (1 + s)* 
= 4f (Fig. 4). 

For points above the parabola, (1 + s)* > 4f, and the 
values of r? and p® are real. We can write Eq. (32) in 
the form 

ae ae : 
> = V(1l + s)* — 4s = Vil + s)}*? — 4 
6¢ 2 ; 
with the understanding that V(1 + s)* — 4s 
V (1 — s)? =1-—s,nots — 1. Itis evident from this 
form that 
1), r? is positiy e 


when s < 1 (below the line s = 


p’ is positive when s < f and negative when s > f. 
For points above s = | we can write Eq. (32) in the 

form 
? - 

6C |: 


= — V(s+ 1)? —4s + V(1+s)? -—4 


with the understanding that V(s + 1)? — 4s 


Vi(s — 1)? =s —1,not1 —s. Itis evident that 


when s > 1, p” is negative 


r? is positive when s > f and negative when s <f. Thus 
by drawing the lines s = 1, s = f (tangent to the para 
bola at 1, 1) the region above the parabola is divided 
into four regions, in each of which the signs of r? and p? 
are known (Fig. 4). 

For points below the parabola, (1 + s)* < 4f, and r* 
and p* are the complex quantities given by 
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Z,ur 
oe ——} 
= x P 
SR comet a 
Fic. 5. Sandwich with free ends 
ms ; z ” 
se =l—-stt1VvV 4 — (1 + s)* (33) 
6c 2 
r 
Then ) = aq, + 1p; 
V 6C ‘op 
where 


m= VVf—-st(1/2)(1-s5, f= 


VVf—s—(1/2)(1—s) (34) 


The characters of ry and p are now known for all 
(positive) values of fand s. The possible combinations 
are 

Region 1: yr and p are both real 

Region 2: and p are complex conjugates 

Region 3: and p are both imaginary 

is real, p is imaginary. 


ss & & 


Region 4: 


(4) THE SANDWICH WITH ‘‘FREE’’ ENDS 


The end conditions considered in this section are indi- 
cated by Fig. 5. The ends of the face sheets are subject 
to horizontal thrusts P(per unit width), and they are 
not displaced vertically (w = 0). No edge shear r, is 
applied to the core—it is excluded by overall equilib- 
rium. The modes considered are of the antisymmetrical 
class of Section 3, with symmetry about the middle 
vertical. 

Since the thrusts P in Fig. 5 are strictly horizontal, 
the shear force Q, which is normal to the deflection 
curve, is given (at each end) by Q = —P dw/dx. In- 
serting this in Eq. (21), with r. = 0 at the ends, and 1/ 
given by Eq. (23), we find 

D(d*w/dx*) + P(dw/dx) = 0 forx = +/ (35) 


The other conditions—no deflection or bending moment 
at the ends—are 


d’w/dx? = 0, forx = +/ (36) 


The mode being symmetrical about the middle vertical 
we take in the general solution Eq. (31) only the terms 
even in x 

w = A»cosh rx + A, cosh px + Ag (37) 
The end conditions, Eqs. (35) and (36), give, in a dif- 
ferent order, 


A. cosh rl + A,cosh pl + Ag = O (38) 
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Aor cosh rl] + Ayp? cosh p/ = 0 39 


A,[r? + (P/D)r] sinh rl + 

A,;(p? + (P/D)p) sinh pl = 0 (40 
The last two yield, as the condition that A», Ay may be 
different from zero 


r?p[p? + (P/D)] sinh pl cosh r/ = 
p’r[r? + (P/D)] sinh 7! cosh pil (4! 


The first possibility here is pp = 0. Then p’r* = 0, and 


from Eq. (32) this means 


fi—s+V(1+s)? — 4f] x 
[’—-s—-V(i+s?-—4f] =0 (42 


which reduces to s = f. But now the differential Eq. 


(28) loses its last term and becomes 
(d§w /dx®) + (12/t?)(S — C)(d*w/dx*) = 0 (43 
so the general symmetrical form for w is now 
w = A,cosh px + Ag + Asx? (44 


where 


2 
-—-—CYf-1) (*& 
p v\ 


The end conditions now yield, as the condition for non- 


vanishing Ay, Ag, As 
p’ cosh pl:(P/D)-1 = [p* + (P/D)p] sinh p/ (46) 


Since P is fixed (s = f), and hence p also [Eq. (45) ], Eq. 
(46) becomes a condition which must be satisfied by 
constants of the sandwich. In general they will not be 
such as to satisfy it, and therefore the possibility of 
satisfying Eq. (41) by pr = 0 can be dismissed. 
Removing the factor pr, Eq. (41) may be written 


r[p? + (P/D)| tanh pl = [r? + (P/D)| tanh r] (47 


To discuss this we consider separately the several re 
gions of Fig. 4. 

Region 1: (r, p, real). 
solve the transcendental Eq. (47) for P, which is im- 
But as / is increased, the hyperbolic 


So long as / is finite we have to 


plicit in p and +. 
tangents approach unity, provided there 7s in the limit 
a finite P which makes p, ry nonzero. The limiting form 
of Eq. (47) is then 
f= [p> + (P Di = plr> + (P D)| or 

[pr — (P/D)|(p — r) = O (AS) 
The possibility p = r requires us to abandon Eq. (31) as 
no longer the general solution. The new solution has 
the symmetrical part 

w = A.cosh rx + A yx sinh rx + A, 

This leads, as did pr = O, to its own critical condition, 
which with p = r fixes a critical load but also requires a 
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relation between the constants of the sandwich. In 
general therefore the buckling will be such that p # r. 
Then Eq. (48) reduces to 


pr = P/D (49) 


The product pr is given in terms of S, C, and F by Eqs. 
29) and (30), these quantities being defined in Eqs. 
15), and so Eq. (49) is easily solved for S. In terms of 


sand f, defined below Eq. (32), the result is 


)+ (1/2) Vi + 4f (50) 
This parabola is shown in Fig. 6, which also reproduces 
the regions of Fig. 4. What has been shown is that the 
part of it within Region 1 gives the critical loads for 
corresponding values of f. To obtain critical loads for 
higher values of f it is necessary to consider Eq. (47) in 
Region 2, which is done in what follows. 
(r, p, complex conjugates). The values of 
Writing 


Region 2: 
r, p are now given by Egs. (33) and (34). 


a = V6C/t a, 8 = V6C/i 6, (51) 
we have 
r,p=axtw@ (52) 


After some reduction the critical condition Eq. (47) is 


Bla? + B? — (P/D)] sinh 2a/ + 
ala? + 6? + (P/D)]| sin 28/ = 0 (583) 


We now seek finite critical loads for large /. Then s, and 


aand @ are finite, and as / increases sinh 2a/—> «. The 
limiting form of Eq. (53) is 

Bla? + B? — (P/D)| = 0 
Vanishing of 8 would mean p = r and is therefore ex- 


The critical condition is thus the vanishing of 
It reduces, by Eqs. (51) and (34), to 


cluded. 
the bracket. 

c+s—f=0 (54) 
This is the same parabola as was found for Region 1, 
and so in Fig. 6 the parabola is continued through Re- 
gion 2. Thus it gives the (lowest) critical loads for all 
values of f, and this completes the analysis for “‘free”’ 
ends. 

Evidently the results are also valid for ‘‘pinned’’ 
ends, since the displacement and bending moment 
vanish at the ends. It is possible to have other “‘pinned- 
ends’ conditions, with nonvanishing vertical end reac- 
Equilibrium then demands shear stress 7, on the 
Any such 


tions. 
ends of the core in the opposing direction. 
atrangement of end forces seems improbable in tests or 


working conditions. 


(5) THe SANDWICH WITH CLAMPED ENDS 
for the face 


If there were 


Fig. 7 indicates the end conditions 
sheets—zero deflecticn and zero slope. 
nonzero shear force Q in the sheets at the ends, it would 


be necessary to assign to 7, on the ends of the core slab 
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Fic. 6. Critical loads for sandwich with free ends 
the value given by 7r.h — 2Q = 0, to ensure overall 
equilibrium. The natural condition to take is that r, 


and QV each vanish. It then follows from Eqs. (21) and 


(23) that d*w/dx* vanishes. The end conditions are 


therefore 


dw/dx = 0, d*w/d*x = 0 atx = +/ 


Taking from the general solution Eq. (31) only terms 
even in x 

w = A,cosh rx + A,cosh px + Ag 
we have from Eqs. (55) the condition that A», Ay, As be 
nonzero as 


(ry? — p*)rp sinh 7/ sinh p/ = O (56) 


The possibilities of 7? — p* or pr vanishing may be dis- 
missed again as requiring special relations between the 
constants of the sandwich. The condition for buckling 
therefore reduces from Eq. (56) to 


sinh 7/ sinh p/ = O Sa) 


For r, p real (Region | of Fig. 4) this can only be satisfied 
by r or p vanishing, one of the cases dismissed above. 
Passing to Region 2 of Fig. 4, and writing 7, p = a + 18 
as in Eqs. (51) and (52), we find that Eq. (57) becomes 


cosh 2a/ — cos 28/ = 0 


and this has no solutions other than a = 8 = O, which 1s 





again within the excluded case pr = 0. 
P P 
A \ 


7 
v0 








P /7 


Fic. 7. Sandwich with clamped ends 
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Considering next Region 3 in which p’, r? are negative, 
we may write 7 = 71y1, p = ty2, and then the buckling 
condition Eq. (57) becomes 


sin y;/ sin yo/ = 0 (58) 


leading to the three possibilities 


(3a): sin y;/ = O, sin yol 0, 
(3b): sin y;/ ¥ 0, sin yo/ = O (59) 
(3c): sin y,/ = Oand sin yl = of 
For (3a) we have y; = mz and therefore since r = 171, 
r? = —n*r’/l*, and from Eq. (32) 
e nx? z : 
— 3 =1l1—s+ Vi(1+s)? — 4f (60) 
oC Tr j 
On squaring out the radical this gives 
tf J Y x 
ase -~2+ §+ where P= 24 
E Ge OF 
(61) 


For finite length / we have to choose the integer ” so as 
to obtain the smallest possible s._ For infinite length we 
can treat //n, and hence é as continuously variable. The 


value of & which minimizes Eq. (61) is 2Vf. The 
corresponding half-wave length //7 is 
rt/¥12C(V f — 1) (62) 


We may therefore accept the minimum of Eq. (61) only 
when f > 1, no new restriction since this inequality pre- 
vails in Region 3. This minimum gives s = 2V f — 1 or 
(s + 1)? = 4f, which is the parabola in Fig. 4 bounding 
Region 2. Here of course we reach it as a limit from 
Region 3. 

Returning to Eqs. (59) and considering the second 
possibility (3b), we find Eq. (60) except for a change of 
sign before the radical. But Eq. (61) is again obtained, 
and we are led to the parabola as before. 

For the third possibility (3c) we have y;/ = uz, y2/ = 
mr with n and m different integers. But then we have 
Eq. (61) as it is, and also with m substituted for x. To 
any given s in Eq. (61) (above its minimum) there are 
two possible values of £, but if the first corresponds to an 
integer m the second in general will not correspond to an 
integer (m). Thus m and m must be the same integer. 
It follows that the deflection w is a pure cosine (Fig. 7). 
This is true also for (3a) and (3b). The sinusoidal as- 
sumption is therefore acceptable here. 

It remains to consider Region 4 of Fig. 4, to find the 
values of s for buckling when f < 1. Here 7 is real, p 


imaginary, and Eq. (57) becomes 
sinh r/ sin yol = O 


Since r = 0 is excluded, we must have sin y2/ = 0 and 
y2ol = nm as under the preceding case (3b). Thus Eq. 
(61) holds again. But now, since f < 1, the minimum 
of Eq. (61) cannot be obtained for-any real half-wave 
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length //n. The lowest admissible value of s from Eq. 
(61) is the lowest we can have without passing out of 
Region 4, and this is s = f—which is higher than the 
minimum of Eq. (61) since f < 1. 

As s approaches f, p approaches zero and so therefore 
does y2, or nw/l. Thus m remains finite as / is in- 
creased—the half-wave length goes to infinity with /. 
Hence for f < 1 we may conclude that the deflection 
curve consists of a single complete cosine wave (crest to 
crest). 

The values of s for buckling with clamped ends (/ in- 
finite), represented by the parabola (1 + s)? = 4f when 


f > 1, and the straight line s = f when f < 1, are shown 


in Fig. 8. This line and curve clearly lie well above the 
parabolic curve for the pinned ends of Fig. 6. The 
symmetrical wrinkling modes of Figs. 2a, b have a com- 


mon critical load S = 2V CF/6 which is s = 2V f/3. 
This parabola lies above the line s = f in Fig. 8 if f< 
4/3. 
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A Theoretical Investigation of the Oscillating 
Control Surface Frequency Response 
Technique of Flight Flutter Testing 


R. A. PEPPING* 
McDonnell Aircraft Corporation 


SUMMARY 


A theoretical background is presented for a frequency response 
technique of flight flutter testing using harmonic deflections of a 
control surface as the excitation mechanism. A graphical method 
is outlined for defining system flutter stability from measured 
driving control surface hinge moment per unit control surface de 
flection. Computed results for a specific flutter system are pre 
sented to exemplify the method of data interpretation. It is 
shown that without complicating the interpretation of test results 
the control surface restraint of the test article can be tailored to 
provide added or artificial stability for testing purposes. This 
indicates the possibility of flight flutter testing, with safety, up 
to or slightly beyond the critical flutter velocity of the normal 
configuration and obtaining measured data from which to define 


this critical velocity 


SYMBOLS 


Coordinates 
h = airfoil elastic bending deflection. See Fig. 1 
a = airfoil elastic torsional deflection. See Fig. 1 
B = control surface rotation. See Fig. 1 
Bo = B, subscript ‘‘O”’ designating out 
B; = series link deflection as equivalent control surface 
rotation. See Fig. 2 
Be = 8; — Bo, displacement error acting on restraining 
impedance. See Fig. 16 
R = reference excitation signal to servo drive. See Fig 
16 
Ha = feedback signal to servo drive 
E = R — Ha, signal error to servo drive 
x; = input displacement to power cylinder valve 
K = output displacement of power cylinder 
\ X; — X,., power cylinder valve deflection 
Control Surface Restraint Terms 
Cg = control surface rotational damping coefficient 
Kg = effective rotational elastic restraint provided by con- 


trol system as seen at control surface, defined by 
Eq. (2.03) 


Ky = effective linear spring rate of control system inboard 
of series link driving mechanism 
K = effective linear spring rate of control system out 


board of series link driving mechanism 


Zg = effective impedance of control system as seen at con- 
trol surface, defined by Eq. (2.02) 

Z, = effective linear impedance of control system inboard 
of series link driving mechanism 

F = power cylinder output force per unit displacement of 


input to valve and zero cylinder motion 


Presented at the Aeroelasticity Session, Twenty-Second 
Annual Meeting, IAS, New York, January 25-29, 1954. 
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Flutter Determinants and Determinant Elements 

F = force or moment partial derivative which, when 
j) + k, is comprised of an aerodynamic and inertia 
term, and when j = & is comprised of an aerody 


namic, inertia, and elastic term 


(Ege Egs — Zg,, control surface hinge moment per unit 
control surface deflection, for zero control sur 
face restraint 

D = stability determinant of the basic system with a 
finite control surface restraint Zg 

D = stability determinant of the basic system with the 
control surface unrestrained 

D = stability determinant of the basic system with the 
control surface restrained by a torsion spring r?K» 

Nees = stability determinant of basic system with infinitely 
rigid restraint of the control surface in rotation 

Nea = minor of the 8a element of the determinant formed 


by the coefficients of the right side of Eq. (3.01 


Transfer Functions 


5 8;/E, driving servo linear displacement per unit 
input signal 

y a«/8, airfoil torsional deflection per unit control sur- 
face rotation 

G a/E, airfoil torsional deflection per unit input sig- 
nal to driving servo 

H = transfer function of stabilizing feedback network 

Miscellaneous 
Ig = control surface moment of inertia about hinge line 
V Kg 

wg = = control surface natural frequency 
138 

w = circular frequency of oscillation 

r = effective length of control surface bellcrank 

g = structural damping coefficient 


(1.0) INTRODUCTION 


ITH THE RAPID FORWARD strides being made in 
\ \ the design of aircraft to achieve high speed and 
the concurrent reduction in the structural stiffness to 
strength ratio of many aircraft components, there has 
been a rapid dwindling of the margin of safety from 
flutter of aireraft designed from strength considerations 
only. These conditions have focused attention on the 
flutter problem, which is no longer a secondary design 
consideration, and emphasized the need for full-scale 
flight flutter testing to preclude unnecessary weight 
penalties resulting from overly conservative theoretical 
flutter considerations. 
Of the various transient and frequency response tech- 
niques which can be employed for full-scale flight flutter 
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testing, one of the most promising from a design and de- 
velopment engineer’s viewpoint appears to be the fre 
quency response technique using an harmonically os- 
cillating control surface to supply the driving forces. 
The frequency response technique provides complete 
excitation of the higher frequency modes sometimes 
not effectively excited by practical transient response 
methods and allows the use of measured data for de- 
fining the system stability which perhaps is more easily 
obtained than the data which would be required with 
the transient response technique. Employing control 
surface motion to drive the dynamic system permits 
using the readily available air as a convenient force 
amplifier and utilizing the airplane control system, 
which is designed for a similar purpose, to transmit the 
driving forces with a minimum of added complication 
and equipment. 

The purpose of this report is to outline the oscillating 
control surface frequency response technique and a 
method of interpreting data measured in flight to define 
the flutter stability of a system. The flight test data 
employed in the method is the driving hinge moment 
per unit control surface deflection required to maintain 
steady-state oscillations. 

The procedure for data interpretation is largely 
graphical and is based on linear theory. The method is 
similar to that of Nyquist’s.* Test results are inter- 
preted so as to provide, when practical, the conven- 
tional flutter stability diagram of velocity versus con- 
trol surface natural frequency (I vs. wg). An alternate 
interpretation provides curves of stabilizing phase mar- 
gin and frequency versus velocity for a specific con- 
figuration of control surface restraint. 

The driving hinge moment per unit control surface de- 
flection is measured at the control surface and is inde- 
pendent of the conditions of control surface restraint 
existing in the test article. It follows, then, that stabi- 
lizing conditions of control surface restraint and con- 
trolled surface motion can be introduced without alter- 
ing the driving hinge moment characteristics. Thus 
artificial stability (or stability for testing purposes) can 
be introduced as a safety measure without interfering 
with the interpretation of the system stability for the 
normal configuration of control surface restraint. 


(2.0) IDEALIZED DyNAMIC SYSTEM 


For the purpose of this theoretical investigation the 
physical system for which the frequency response to 
control surface harmonic oscillation is to be obtained is 
idealized as an airfoil having three elastically uncoupled 
degrees of freedom consisting of airfoil bending and 
torsion and control surface rotation. See Fig. 1. The 
addition of other degrees of freedom would complicate 
the analysis without adding to the basic concepts to be 
derived from the study. 

The control surface of the system shown in Fig. | is 
considered to be subjected to a driving hinge moment 
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Pigure 1. - Idealized Degrees of Freedom 
of Flutter System 
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Figure 2. - Idealization of Control System 


A; which ts applied through the airplane control system 
by a series link drive. Without considering the details 
of the driving mechanism, the combination of control 
surface, control surface restraint, and driving displace- 
ment or hinge moment can be idealized in the manner 
indicated by Fig. 2. 

The equation relating the restraint of the control sur- 
face and the applied hinge moment can be written 


Z3(8: — 8) = M, (2.01 


where Z, is the effective impedance of the control system 
as seen at the control surface and is defined by the ex- 
pression 


(2.02 


Zp = 2,Kor?/Z, + Ko 


The quantity Z, is the effective impedance of the con- 
trol system inboard of the displacement device and may 
be a function of the frequency of oscillation. When Z, is 
a simple spring A,, the impedance of the control system 
is a torsion spring A, defined by the expression 


Kg = K,Kor?/(K, + Ko) 


(2.03 


From Eq. (2.01) it follows that when the input con- 
trol surface driving deflection 6; vanishes, the external 
hinge moment — ./, is resisted by the internal restrain- 


ing moment 7-8. 


(3.0) EQUATIONS OF MOTION 


Employing linear theory the oscillatory motion of the 
idealized dynamic system can be prescribed by the 
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three simultaneous equations, 
Ws = (Egs)0°8 + Eg,:h + Ee,- a (a) 
O = Eng B+ Em h + Eng’ (b) 
FEap°B ss Eanth + Fan’ & (Cc) 


(3.01) 


where /,, is a force or moment partial derivative which, 
when j} ¥ k, is comprised of an aerodynamic and inertia 
term and, when j = k, is comprised of an aerodynamic, 
inertia and elastic term. 

Since the applied hinge moment .\/, is transmitted 
through the control system which would otherwise 
supply the control surface restraint, the derivative 
Ez3)) does not include a restraint term, i.e., 

(Egs)0 = Egg — Ze (3.02) 

Solving the set of simultaneous equations, Eq. 

3.01), for 3/8 by the method of determinants yields 


the expression 


M,/B = Do/Nees (3.03) 
where 
D) = the stability determinant of the basic system 
with the control surface unrestrained (free 
to rotate) formed by the coefficients of the 
right side of Eq. (3.01) 
Naz = the minor of the 88 element of the determi- 


nant of coefficients; but more significant is 
the stability determinant of the basic sys- 
tem with infinitely rigid restraint of the 
control surface in rotation. 

From Eqs. (3.01) and 3.02) it can be seen that 


D = Dy + Zg-Nas (3.04) 


where D is the stability determinant for the system with 
the control surface restrained by Z3. 

Adding Z, to both sides of Eq. (3.03), reducing the 
right side of the expression to a common denominator 


and substituting Eq. (3.04) yields 
7 


(M,/8) + Zs, = D/Nee 


or (3.05) 
(M,/Z,8) + 1 = D/ZgNee ual 
By a similar manipulation, 
(B/M,) + (1/Z,) = D/ZgDo | 
or (3.06) 
(Z,8/M,) + 1 = D/Do ol 


Both Eqs. (3.05) and (3.06) define the same system. 
The control surface restrained stability is shown to be 
related to the measurable driving hinge moment per 
unit control surface deflection which is required to 
maintain steady-state oscillation. As will be discussed 
later, there are advantages to each expression in the in- 
terpretation of test results, depending on the character- 
istics of the system for the test conditions. 
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(4.0) INTERPRETATION OF TEST RESULTS 


4.1) Velocity Versus Control Surface Natural Frequency 
V vs. ws) Diagram 


(4.1.1) Mg/8 Plot.—At a constant value of the for 
ward velocity |’ the frequency w of the driving oscilla- 
tory hinge moment can be varied and the vector 1/,/8 
determined from flight test measured values of applied 
hinge moment and control surface deflection. If the 
test velocity is sufficiently high, the curve formed by the 
tips of the vectors ./,/8 when plotted on the complex 
plane might appear approximately as shown in Fig. 3. 

It is seen from examination of Fig. 3 that when a real 
vector of value Ag, or Ag, is added to the 1/7, 8 vectors, 
the curve formed by the tip of the combined vector, 
which from Eq. (3.05) is 


D/Nog = (Mg/B) + Kg 


passes through the origin. This is equivalent to the sta- 
bility determinant D vanishing (Vg, would remain 
finite), which means that the system would be neu- 
trally stable for these two values of control surface elas 
tic restraint. 
neutral stability occurs can also be obtained from the 


The frequency of oscillation at which 


plot. 

For an intermediate value of Ay, between A,, and 
Kg,, there remains a negative imaginary component of 
the combined vector when the real portion is zero. 
The curve formed by the nose of the vector (17,/8 + Kg) 
would pass the origin in the manner indicated by Fig. 4. 
In terms of work being done, the unbalanced negative 
imaginary vector indicates that for this value of Ay, 
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Figure 5. - Constant Velocity, Varying Frequency 
Plot of Vector Mg/p 
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Plot of Vector (Ne/p + Kp 
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Figure 5. - Velocity Versus Control Surface Natural 
Frequency Stability Diagram 











energy must be absorbed to maintain steady-state 
oscillation or the system would be unstable. It follows 
then that for the configuration being considered insta- 
bility exists for the conditions A,, << A, < Ag,. The 
two values of A, thus defined which yield neutral sta- 
bility, plus the knowledge that intermediate values 
produce instability, form part of a conventional velocity 
versus control surface natural frequency diagram. By 
testing at successively higher velocities, and repeating 
the analysis of data, the complete diagram can be ob- 
tained experimentally and might appear approximately 
as shown in Fig. 5. The stability of the specific con- 
figuration being tested is, of course, defined by the 
actual value of w, for the system. 

In the discussion so far, the tacit assumption has been 
made that at the speeds being tested the control surface 
fixed stability of the system is positive, 1.e., all the roots 
of Ngg are stable or the test speed is below the asymp- 
tote 1’ = A shown in Fig. 5. Under this condition the 
analysis of test results is not clouded by the apparent 
discontinuity which would be produced by Ng, vanishing 
and the following general rules for data interpretation 
can be formulated. 

Rules for interpreting 1/,/8 Plot (all roots of N43, 
stable) : 

1. From measured data make a complex plot of 
M,/B. 

2. Form the M/,/8 + Kg plot by adding a real vec- 
tor A, to the varying complex /,/8 vector -a positive 
real vector if the assumed restraint is positive and a 
negative real vector if the assumed restraint is negative. 

3. If the curve formed by the tip of the 17,/8 + Kz 
vector passes through the origin, the system would be 
neutrally stable if the control surface were restrained by 
this value of Kg. 

4. If the origin is between successive crossings of the 
real axis by the curve formed by the 1/,/8 + Kg vector 
and the imaginary component of 1/,/8 + Kg is negative 
when the real portion is zero, the system would be un- 
stable if the control surface were restrained by this 
value of Ky. 


-AUGUST, 1954 

There will be occasions when it will be desirable to de- 
fine the stability boundary over a speed range in which 
the control surface fixed stability is not positive. In 
Fig. 5, this would be the speed range from A to B. Re- 
ferring to Fig. 3, it is evident from the nature of the 
M,/8 plot that the second crossover of the real axis will 
approach minus infinity as the critical speed for control 
surface fixed (_V = A in Fig. 5) This 
occurs as the denominator of the right side of Eq. 
(3.05(a)), gg, approaches zero. 
when V approaches A, designated (—Kz,,, j0), will 
occur with the vector ./,/8 rotating clockwise as shown 
in Fig. 6. Above the critical velocity, V = A, the 
second crossover on the negative real axis does not 


is approached. 


The second crossover 


occur; however, the vector 1/;/8 may cross the positive 
real axis while sweeping counterclockwise. See Fig. 7. 

The two crossovers, one on the negative real axis and 
one on the positive real axis, are valid indications of con- 
ditions for neutral stability and, having followed the 
trend with increasing speed, it is apparent that the 
complete system would be unstable for aa elastic re- 
straint greater than Ay, out to infinity, and back from 
minus infinity to A,,. The simple rules for interpreta- 
tion of test results previously formulated, however, do 
not in this case provide a clear picture of the conditions 
for stability because of the crossing of the negative 
imaginary axis by the 1/,/8 vector. This crossing does 
not indicate instability with zero elastic restraint, but 
demonstrates that the stabilizing influence of the con- 
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Figure 8. - Constant Velocity, Varying Frequency 
Plot of Vector B/Mg 





trol surface motion for zero elastic restraint can be ef- 
fectively eliminated by sufficient restraint provided 
through a pure damper. 

(4.1.2) B/M, Plot.—Because of the large values of 
M,,/8 and possible confusion in interpreting the 7, 6 


plot near = A, it may prove expedient to employ a 
8/\M, plot. With this diagram the difficulties around 
Nag = O or V = A are circumvented; however, for the 


8/\l, plot to be a practical alternate to the 7,/8 plot 
the control surface free stability defined by Do must be 
That Dy has no unstable roots can be deter- 
It is expected that in 


positive. 
mined from the 1/3/68 diagram. 
the practical case the control surface free stability will 
be positive for speeds above V = A, if it is of interest to 
investigate this speed range. The plot of 8 \/, cor- 
responding to the plot of Fig. 7 would appear approxi- 
mately as shown in Fig. 8. It can be seen from ex- 
amination of the plot that when a real vector of value 
1/Kz,, or 1 Ag, is added to the 6/17, vector, the com- 
bined vector, which from Eq. (3.06) is 


D/KgDo = (8/Mg) + (1/Kg) 


passes through the origin. This is equivalent to the 
stability determinant D vanishing (A,4)) would remain 


finite), which means that the system would be neutrally 
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Figure 9. - Constant Velocity, Varying Frequency 
Plot of Vector 6/Mg + 1 
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Figure 10. - Constant Velocity, Varying Frequency Pict 
of Vector Mp/p, Showing Phase 
Margin at |Mp/pl = Ks 
stable for these two values of control surface elastic 
restraint. 

For an intermediate value of 1/A,, between 1, A, 
and 1/K,, and a frequency above that at which the 
first crossover occurs, there remains a positive imagi- 
nary component of the combined vector when the real 
portion is zero. The curve formed by the nose of the 
vector 8/11, + 1/Ag would then pass the origin in the 
manner indicated by Fig. 9. The unbalanced positive 
imaginary component of 8/ \/, is effectively a negative 
imaginary component of .\/,/6 and, again, from the work 
being done to maintain steady-state oscillation indicates 
that the system would be unstable for this value of 


1/Ke. 
] Kg, > ] Kz > 1/Kg, or 


Instability exists then for 
Ky, < Kg < Kg, 


This information is similar to that obtained from the 
M,/8 plot and is applicable to the determination of the 
V versus wg diagram. It should be noted that the un- 
stable root of Ng, is indicated by a clockwise encircling 
of the origin on the 8/17, plot. 

The rules for interpreting the 6 .\/, plot are the same 
as the rules for interpreting the .\/, 8 plot with the inter- 
change of 8/M, for W,/8 and 1/Kg for Kg, except that 
instability is indicated by an unbalanced positive 
imaginary component of 8 1/3; + 1/A, when the real 
portion is zero. 

Sample computed curves of 17,/8 for a specific prob- 
lem are presented as Appendix A. 

4.2) Margin of Stability 

The interpretation of test results to define a !° versus 
wg diagram will not be, for all test configurations, an 
effective method of determining system stability. In 
some instances it will be necessary to follow the trends 
in flutter stability with increasing velocity over a speed 
range in which the test configuration is stable for all 
values of control surface elastic restraint or of a con- 
figuration for which the control surface restraint cannot 
be represented as a simple spring. Several alternates 
to the previously discussed method of data interpreta- 
tion are outlined for use under these circumstances. 
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Figure ll. - Plots of Stabilizing Phase Margin 
and Frequency versus Velocity 


(4.2.1) Simple Spring Restraint.— For the case where 
the complex /,/8 vector has a positive imaginary com- 
ponent throughout the frequency sweep, the plot of this 
vector on the complex plane might appear approxi- 
mately as shown in Fig. 10. If the control surface re- 
straint can be idealized as a simple spring Ag, the point 
along the 73/8 or 6/M, plot at which 


|\M,/8) = K, or  |B/M,| = 1/Kg 


defines the condition wherein the external and internal 
forces are in balance except for a phase angle dis- 
crepancy. The minimum value of this phase angle 
which is between the /,/8 or 8/M, vector and the 
negative real axis can be employed as a measure of the 
system stability, providing V,, has no unstable root if 
the angle is defined on a 4/,/8 plot or Dy has no unstable 
root if the angle is defined on a 8/M, plot. Using feed- 
back control terminology, an angle measured from the 
negative real axis to the 1/,/8 or 8/ MJ, vector will be 
called a phase margin which, on the 1/,/8 plot is con- 
sidered positive when measured clockwise, and on the 
8/M, plot is considered positive when measured 
counterclockwise. The minimum phase margin for the 
particular condition wherein the moduli of external and 
internal forces are equal will be referred to as the 
stabilizing phase margin and designated 6; see Fig. 10 
for example of graphical interpretation. Subject to the 
previously stated restrictions concerning an unstable 
root of Ng, or Dy: systems having positive stabilizing 
phase margin are stable, whereas systems with negative 
stabilizing phase margin are unstable. 

At each test velocity the stabilizing phase margin and 
the associated frequency w can be determined. Plots of 
6 and w versus V, which might appear approximately as 
shown in Fig. 11, can be constructed and used to effec- 
tively follow the trends in stability with increasing 
velocity. 

(4.2.2) Power Cylinder Control Restraint.—In some 
test configurations an hydraulic power cylinder will be 
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employed in the control system somewhat in the manner 
indicated by Fig. 12. 

The output displacement of the power cylinder X, is 
For the 
system shown, the valve displacement is defined by the 


controlled by feedback to the hydraulic valve. 
expression 
X, = 
or (4.0] 
X;, = X,+ Xo es 


Though not necessarily true for all installations, the in- 


put displacement to the valve X; is considered to be 
unaffected by structural deformation and motion of the 
airplane. 

The characteristics of a power cylinder hydraulic 
valve are known to be nonlinear; however, reasonable 
success has been attained in analyzing the stability of 
power control systems using linear theory and defining 
the stability for small perturbations about the various 
quasi-steady-state operating conditions which are phys- 
ically realized.° This technique is considered appli- 
cable to the interpretation of the flight test measured 
data. 

The linear equation of the power cylinder total out- 
put force F can be expressed in the form 


f= FX; = Bike (4.02) 


Both of the partial derivatives /; and —Z, are charac- 
teristics of the power cylinder plus power cylinder 
backup structure and are functions of the frequency of 
oscillation. The linear idealization of the system of 
Fig. 12 can be represented in the manner shown by Fig. 
is. 

The balance of the oscillatory hinge moments for the 
system of Fig. 13 can be defined by the two equations, 


Mz + 2,8 = F(X, + X,) (a) 


Ky 
K.+2Z, (4.03) 


Mg + r°K2oB = rK2X, (p)) 
where Z, is the effective impedance of the control system 
as seen at the control surface and is defined by Eq. 
(2.02). The stability of the system is determined by 
the left side of Eq. (4.03(a)). If the vector sum 1/3 + 
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Pigure 12. - Power Cylinder - Control 
Surface Combination 
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Effective impedance 
of power cylinder 
plus power cylinder 
backup structure. 





Fy -% 


Cylinder output force 
produced by input dis- 
placement of hydraulic 
valve. 


Figure 13. - Linear Idealization of Power 
Control System 
Z,°8 is zero, the system will oscillate without the as- 
sistance of displacement of the input rod to the hy- 
draulic valve. 

Considering only the left side of Eq. (4.03(a)) the re- 
lationship between the measurable .1/,/8, control system 
impedance and the system stability is defined by ex- 
pression Eqs. (3.05(b)) or (3.06(b)). Since Z, is con- 
sidered to be a complex vector which varies with fre- 
quency, this form of the equation will be more con- 
venient for graphical interpretation than Eqs. (3.05(a)) 
or (3.06(a)). 

Having measured the vector .1/;/6 for a constant test 
velocity and varying frequency of control surface os- 
cillation, and having measured or computed the con- 
trol system impedance Z, for a particular steady-state 
operating condition, the product vector J/,/2,;0 or 
Z;3/ M, can be computed and plotted on the complex 
plane. The condition of neutral stability is defined on 
the plot by the product vector passing through the 
point (—1, 70) on the real axis. When this occurs the 
right sides of Eqs. (3.05(b)) and (3.06(b)) vanish which 
is equivalent to the stability determinant D vanishing 
(Z3, Ngg, or Do remain finite). The plots of the product 
vectors are interpreted in the same manner as outlined 
in Section 4.2.1 for the respective 1/,/8 or B/ \/, plots, 
except that the strategic point on the real axis is 
(—1, 70) instead of (—Kg, 70) or (—1/Kg, 70). The 
stability of the system is indicated by the phase margin 
at the point where the product vector crosses a ref- 
erence circle of unit radius with its center at the origin, 
providing V,, has no unstable root if the phase margin 
is defined on a 73/238 plot or Dy has no unstable root 
if the angle is defined on a Z7,3/.\/3 plot. The presence 
of an unstable root of N33 or Dy can be determined from 
the plot of the 8/\/, or \73/8 vector. 

Employing the 1/,/Z,8 and/or the 738, \/, plot, the 
stabilizing phase margin 6 and frequency w at which the 
unit circle is crossed can be determined graphically from 
measured variations of the vector 1/,/8 obtained at 
successively higher test velocities and the known linear 
characteristics of the control system impedance Z, for 
various applicable steady-state operating conditions. 
Plots of 6 and w versus velocity, similar to those of Fig. 
11, can be constructed and used to follow the trends in 
stability with increasing velocity. 
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When the oscillatory deflection of the control surface 
is accomplished by a displacement device located be- 
tween the power cylinder and the control surface, the 
effective impedance of the power cylinder plus backup 
structure can be measured in flight during a frequency 
response test. With the input displacement to the 
valve X; held constant, the ratio of the oscillatory load 
at the output end of the power cylinder to the oscilla- 
tory output displacement of the power cylinder can be 
measured. This ratio, within the limits of the assump- 
tion of linearity, is the effective impedance Z,. 

When control surface motion is accomplished by har- 
monic displacement of the input rod to the power cylin- 
der, the trend in system stability can also be determined 
from the measured ratio of valve displacement to power 
cylinder output displacement. Here, again, the con- 
cept is dependent on the applicability of linear theory. 
The relation between XX, and system stability can be 
obtained from Eqs. (4.03(a)) and (4.03(b)) and takes 
the form 
X, K,+ 2, 


l = 
pe ™ Fy 


(13/8) + Zs 
(.7,/8) + r°Khe 


(4.04) 


Making the substitutions indicated by Eq. (3.05(a)), 


Eq. (4.04) becomes 


X,/Xo +1 = [((K2 + 2:)/Fi]:(D/D2) (4.05) 
where 
D = system stability determinant 
D. = stability determinant of the system with the 


control surface restrained by a torsion spring 

r°Ko. 
From measured values of XY, X, 
stant test velocity, a plot of Y,/X, 
The condition of neutral 


obtained at a con- 
on the complex 


plane can be constructed. 


stability is defined on the plot by the vector XY, X, pass- 
ing through the point (— 1, 70) on the real axis. Subject 


to the restriction that D» has no unstable root at the test 
velocity, the phase discrepancy between X , and Y, when 
X,/Xo| = 
system stability. 
D. can be determined from a plot of the 1/3 8 vector. 


1 can be employed as a measure of the 
The presence of an unstable root of 


(5.0) SAFETY DEVICES 


Since the quantity ./,/8 which is used for stability 
determination is measured at the control surface, the 
introduction of any devices inboard of this point in the 
control system will not alter 1/,/8. Utilizing this fea- 
ture, the restraint of the control surface for 8; = 0, the 
driving mechanism at rest, can be tailored through use 
of an hydraulic damper and/or a stabilizing feedback 
circuit, in addition to a favorable adjustment of the 
elastic restraint, so as to provide added stability for 
testing purposes without interfering with or complicat- 
ing the interpretation of test results. For the initial 
phase of a test the control surface restraint can be 
tailored on the basis of a theoretical analysis and as the 
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Figure 14. - Idealization of Control System 
with Hydraulic Damper Installed 
as a Safety Device 

test progresses to higher velocities this restraint can be 
modified in accordance with the changes indicated by 
analysis of measured test results. By this procedure, 
flight testing can be carried to speeds beyond the critical 
velocity of the airplane in its normal configuration while 
maintaining a positive margin of stability against 
flutter. The use of an hydraulic damper and the use of 
a feedback circuit as a safety device is discussed in the 
following section. 

5.1) Hydraulic Damper Installation 


For the particular case where the flight flutter testing 
is to be confined to velocities for which the control sur- 
face fixed stability is positive, i.e., velocities for which 
N,; has no unstable roots, probably the simplest safety 
device which can be employed is the installation of an 
hydraulic damper in the control system. The damper 
should be introduced as a parallel link just inboard of 
the point where the control surface driving hinge 
moment is measured, as indicated by Fig. 14. 

The equation relating the restraint of the control 
surface and the driving hinge moment can be written 


Ky: (8: — B) — jCwB = M, (5.01) 


Since the effective impedance of the control system 
for B; = O, the driving mechanism at rest, is Kz + jC, 
it follows from Eq. (3.05(a)) that 


D/Npg = (M3/B) + jCgw + Ke (5.02) 


where D is the stability determinant for the system with 
the control surface restrained by the spring plus damper, 
Kg + Cg. 





A Imaginary 











Figure 15. - Constant Velocity, by j 
Frequency Plot of Vector (Me/p + JCpa, 
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For a fixed value of A, the characteristics of the vee- 
tor 1/33 + jCgw define the stability of the system with 
damper installed. 

The combined vector J/,/8 + j7C3w can be measured 
by instrumentation located just inboard of the hy- 
draulic damper while the vector 1/;/8 can be measured 
on the same test article by instrumentation located out- 
board of the hydraulic damper. Plotted on the complex 
plane, the measured variation in these two vectors with 
increasing frequency measured at a constant test ve- 
locity, might appear approximately as shown in Fig. 15, 

Employing the methods previously outlined, the 
measured values of 17,/8 + jCgw and \/,/8 which are 
obtained at successively higher test velocities can be 
interpreted to define the trends in stability with in- 
creasing velocity for both the system with and the sys- 
tem without hydraulic damper. If the critical mode 
being investigated involves considerable control surface 
motion the system with damper installed would be ex- 
pected to have the higher critical velocity. Under these 
conditions, flight flutter testing could be conducted, 
with safety, up to or slightly beyond the critical flutter 
velocity of the normal configuration without damper in 
order to obtain measured data from which to define this 


critical velocity. 
5.2) Stabilizing Feedback Circuit 


When it is desirable to test at speeds approaching the 
critical flutter velocity for the fixed control surface con- 
figuration the hydraulic damper may prove ineffective 
as a Stabilizing device. Under these circumstances, 
however, controlled motion of the control surface could 
be used to produce added or artificial stability in the 
test article. 

If, in the test article, electrical pickups are provided 
to measure motion of the primary surface, say in the 
torsion degree of freedom, and signals from these 
pickups are fed back to the excitation mechanism so as 
to control the rotation of the control surface in ampli- 
tude and phasing, motions of the control surface could 
be produced which would effectively stabilize the dy- 
namic system. This method of using controlled surface 
motion for stabilization purposes has been employed 
effectively in the field of rigid airplane dynamic stabil- 
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Figure 16. - Block Diagram of Stabilizing Feedback Circuit, Using 
Motion in the Torsion Degree of Freedom to Provide the 
Feeavack Signal to the Control Surface Actustor 
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OSCILLATING CONTROL SURFACE 


ity... Fig. 16 shows roughly in block diagram form the 
feedback circuit which would be required. The driving 
mechanism is considered to be an electrical-mechanical 
or electrical-hydraulic servo which converts a mixed 
electrical signal, consisting of the driving and stabilizing 
feedback components, into mechanical displacement for 
control surface actuation. When the driving signal is 
zero the control surface is effectively geared to the 
primary surface to stabilize the combination. One 
difficulty with this system will be, of course, obtaining a 
sufficiently fast servo that has a flat response throughout 
the test frequency range and can supply sufficiently 
large control surface deflections. 

The inner loop of the block diagram of Fig. 16 repre- 
sents the system of Fig. 2. The subscript zero has been 
added to 8 to indicate ‘‘out,’’ however, 8) = 8 as pre- 
viously employed. The transfer functions for the 
servo, airplane, and feedback circuit, designated S, Y, 
and //, respectively, are the output to input ratios for 
these components of the complete loop. The transfer 
function Y for the system of Fig. 1 can be defined from 
Eq. (3.01) by the expression 


VY = a/B = N,/N (5.03) 


“"oe BS 


where V,, is the minor of the Sa element of the deter- 
minant formed by the coefficients of the right side of 
Eq. (3.01). 

The block diagram for the outer loop of the system of 
Fig. 16 is shown in Fig. 17. The transfer function desig- 
nated G can be defined from the component transfer 
functions by the expression 


G = SY/[(M,/Z,8) + 1] (5.04) 


Substituting Eqs. (3.05(b)) and (5.03) in Eq. (5.04), the 
inverse transfer function is defined by the expression 


1/G = (D/Z,N,,) 1/S (5.05) 


where D is the stability determinant of the system 
without feedback. 

The stability of the complete system of Fig. 16 is the 
stability of the outer loop and can be determined from 
the steady-state input to output ratio R/a. From the 
block diagram of Fig. 17 it can be seen that 


R/a = (1+ GH)/G = (1/G) +H (5.06) 


But 1/G is defined by Eq. (5.05) which, when substi- 
tuted in Eq. (5.06), results in the equation 


R/a = [D + (HSZ,Nog,)]/SZNa, (5.07) 


The numerator of the right side of Eq. (5.07) is the 
stability determinant for the system with the a degree 
of freedom controlled by a feedback circuit to the con- 
trol surface actuator which differs from the stability 
determinant of the system without feedback by the 
factor HSZ,N,.- 

The variation with frequency of the ratio of the ref- 
erence excitation or input signal to the output of the 
controlled degree of freedom (R/a for the case con- 
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Pigure 17. - Block Diagram of Outer Loop 
of Stabilizing Feedback 
Circuit 


sidered) can be measured at each test velocity and, 
applying servomechanism theory, the complex plot of 
this closed loop inverse frequency response can be 
analyzed to define the stability of the controlled system 
and its variation with increasing velocity. The addi- 
tion of the feedback circuit does not alter the character- 
istics of the inner loop of Fig. 16 and the stability of the 
system without feedback can be determined from 
measured values of 1/,/8. 

It is expected that by tailoring the feedback circuit 
to the particular flutter configuration to be flight tested, 
the critical flutter velocity of the test article can be 
made higher than the critical velocity of the normal con- 
figuration. Under these conditions, flight flutter testing 
could be conducted, with safety, up to or slightly be- 
yond the critical flutter velocity of the normal con- 
figuration without feedback in order to obtain measured 
data from which to define this critical velocity. 


(6.0) CONCLUDING REMARKS 


A theoretical background is presented for a frequency 
response technique of flight flutter testing using har- 
monic deflection of a control surface as the excitation 
mechanism. It is shown, on the basis of linear theory, 
that the flutter stability of a system can be determined 
from the oscillatory hinge moment per unit control sur- 
face deflection required to maintain steady-state oscilla- 
tions, which can be readily measured in flight. A 
graphical method is outlined for interpreting the 
measured hinge moment parameter to define system 
stabilty. It is also shown that with this technique of 
testing it is possible to tailor the control surface re- 
straint of the test article to raise the critical flutter ve- 
locity of the test article above that of the normal con- 
figuration without complicating the interpretation of 
test results. With the test article stabilized by an hy- 
draulic damper in the control system and/or a feedback 
circuit which provides controlled motion of the control 
surface, data can be measured in flight from which the 
stability of the normal configuration as well as the 
stability of the test article can be determined by direct 
graphical interpretation. Providing the.stabilizing de- 
vice is effective, it is possible to conduct flight flutter 
tests, with safety, up to or slightly beyond the critical 
flutter velocity of the normal configuration and obtain 
measured data from which to define this critical velocity. 
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The graphical method which is presented for inter- 
preting test results also provides a useful tool for 
theoretical analysis when it is desirable to tailor the 
control surface restraint to an optimum configuration. 
The method is particularly convenient for evaluating 
the benefits of a complex restraint which might result 












































































































































from the use of a spring-damper combination or hy- 
draulic power controls. 
RUDDER ROTATION - FIN TORSION - 
COUPLED FIN BENDING AND AFT FUSELAGE TORSION 
FOR g = 0 ALT. = 10,000 FT. 
sts Pa 
g = —— 
& ‘YY ATL 
; N | 
& ri 
< 
k 400 
2 
a is 
z 
' 200 
0 
0 40 80 120 160 200 2ko 280 320 
@ - RAD. /SEC 
Figure A-1 
Critical Flutter Speed Vs. 
Rudder Rotational Frequency 
—— V._ = XO knots 
—am= Ye © WO hunts 
ee 
30000 — 
hi 7 
\ ' o/ 
' 
j Ps 
20000 @ yy ' / ye fou 
YW ! / f 
1. ZEMEL I db 
~ 10000 xt ° 
: ' ae a 
: ae a 
& es Sp ew oe - as 
: ° —+-— ~ — en 
: a Fy 
i Se « 
= ’ 
‘S  =—«--10000 va 
gz me 
? 
-20000 
-40000-« -30000 -20000-—=s_ 10000 ° 10000 20000 30000 
0% /e)peay - tt-1b./rad. 
Figure A-2 
Plot of Computed Mg /g 
° ¥, = 5CO knots 
Vv V_ = 700 knots 
O v, = 900 knots 
8000 \ T 
4000 ~\ 
: \ ‘y 
2 ‘ a te 
é 
3 es = 
i 4000 Se a 
- 8000 
of 
VA 
or -24000 -20000 -16000s -12000 8000-4000 ° 4000 


(Mg /p)Read - ft-1d./red. 


Figure A-3 
Plot of Computed Mg /p 


AERONAUTICAL 


SCIENCES AUGUST, 1954 
REFERENCES 


! Theodorsen, Theodore, General Theory of Aerodynamic Inst 
bility and the Mechanism of Flutter, NACA TR 496, 1935 


2 Nyquist, H., Regeneration Theory, Bell System Techn. Jour., 


Vol. XI, No. 1, January, 1932. 

3 MacColl, LeRoy A., Fundamental Theory of 
nisms; D. Van Nostrand Company, Inc., New York, 1945 

* Chestnut, Harold, and Mayer, Robert W., Servomechanism 
and Regulating System Design; Vol 
New York, 1951. 

5 Report of the First Piloted Aircraft Powered Surface Control 
System Symposium, Bureau of Aeronautics, Airborne Equipment 


Division, Department of the Navy, Washington, D.C., October, 


1949. 

6 Timoshenko, S., Vibration Problems in Engineering; D 
Nostrand Company, Inc., New York, 1937. 

7 Den Hartog, J. P., Mechanical Vibrations; McGraw-Hill Book 


Van 


Company, Inc., 1940 

8’ White, Roland J., Investigation of Lateral Dynamic Stability in 
the X B-47 Airplane, Journal of the Aeronautical Sciences, Vol. 17, 
No. 3, pp. 1383-148, March, 1950 


APPENDIX 


PLOTS FOR A THREE-DEGREE-OF- 


FREEDOM SYSTEM 


CoMPUTED M,/6 


To provide a picture of some of the variations in the 
M,/8 plot which can be expected, the driving hinge 
moment per unit control surface deflection has been 
computed for a three-degree-of-freedom system using 
incompressible aerodynamic derivatives. The condi- 
ditions for neutral stability are shown to be the same as 
predicted by the I’ vs. g method of solution. 

From successive | vs. g solutions of the equations of 
motion for the three-degree-of-freedom system com- 
prised of fin torsion, coupled fin bending and aft fuselage 
torsion, and rudder rotation for various assumed values 
of rudder elastic restraint, the V vs. ws, diagram shown 
as Fig. A-1 was obtained. For three conditions of for- 
ward velocity, V, = 500, 700, and 900 knots, the driving 
hinge moment for increasing frequency was computed 
for the same dynamic system. These computed results 
are presented in Fig. A-2. Examining the plot, the in- 
fluence of the approaching Vz, solution is apparent at 
the higher frequency end of the 500- and 700-knot curves 
as a large clockwise loop. For the curve corresponding 
to 900 knots, which is above the critical velocity for 
control surface fixed, the influence of the unstable solu- 
tion of NVgz is shown by the large loop formed by the 
M,/8 vector sweeping counterclockwise. 

The section of Fig. A-2 in close proximity to the 
origin has been enlarged and is presented as Fig. A-3. 
From this diagram the values of A, for neutral stability 
can be determined and from these values the correspond- 
ing control surface natural frequencies can be com- 
The moment of inertia for the rudder is 0.415 
The values of w, defined in this manner are 


puted. 
slug (feet). 
in complete 
the V vs. g 
A-l. 


agreement with the values obtained from 
method of analysis and presented in Fig. 
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On the Concept of Stability of Inelastic 
Systems 


ey 


D. C. DRUCKER* 


AND E. T. ONAT?t 


Brown Unwersity 


ABSTRACT 


Simple models are employed to bring out the large and im- 
portant differences between buckling in the plastic range and 
classical elastic instability. Static and kinetic criteria are com- 
pared and their interrelation discussed. Nonlinear behavior in 
particular is often found to be the key to the physically valid 
solution. The nonconservative nature of plastic deformation in 
itself or in combination with the nonlinearity requires concepts 
not found in classical approaches. Conversely, the classical 
linearized condition of neutral equilibrium is really not relevant 
in inelastic buckling. Plastic buckling loads are not uniquely 
defined but cover a range of values and are often more properly 
thought of as maximum loads for some reasonable initial imper- 
fection in geometry or dynamic disturbance 

The models indicate that basically the same information is 
obtained from essentially static systems by assuming initial im- 
perfection in geometric forms as by assuming dynamic disturb- 
ances. One approach complements the other and both are help- 
ful in obtaining an understanding of the physical phenomena. 


INTRODUCTION 


[* THE ANALYSIS OF structures, as in most branches of 
engineering and science, problems are solved by 
simplifying them enormously. The geometry of the 
structure, the material of which it is composed, and 
the loads applied are all strongly idealized. Justifi- 
cation of any idealization requires an analysis of the 
real system and an interpretation of the results. Should 
reasonable deviations produce large changes in the 
result obtained from the idealized system the idealiza- 
tion is not permissible. 

For example, in the basic elastic buckling problem of 
the Euler strut, the assumptions made are that the 
column is perfectly straight, that the load is static and 
is applied along the centerline, and that the material 
is homogeneous, linearly elastic, and free of initial stress. 
None of these are strictly true. 

The usefulness of the Euler critical load computation 
lies in the fact that the calculations based on real col- 
umns show the critical load to be a reasonable limiting 
load as long as the stresses induced are below the elastic 
limit and the imperfections are small, Fig. 1. 


Presented at the Structures Session, Twenty-Second Annual 
Meeting, IAS, New York, January 25-29, 1954. 

* Chairman, Division of Engineering. 
B. Jewett Fellow in Applied Mathematics. 

t An interesting discussion of this general question is contained 
in N. J. Hoff’s Wright Memorial Lecture, Buckling and Stability, 
Journal of the Royal Aeronautical Society, Vol. 58, No. 517, 
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Similar types of calculations have proved equally 
successful for other structural elements where small 
the not 
produce large differences in behavior. 
successful analysis of linearized, idealized, and simpli- 
It seems to have led 


deviations from idealized condition also do 


In a way, this 


fied systems was unfortunate. 
to the idea that the classical linearized theory with its 
static criterion of neutral equilibrium solved the real 
problem. As an actual imperfect system is so trouble- 
some, even in the elastic range, it is easy to see how 
such a situation could develop. When low values were 
obtained experimentally for the buckling loads of 
cylindrical and spherical shells the warning was not 
taken seriously in general. Donnell did attribute the 
discrepancy to initial imperfection but von Karman and 
Tsien' explained the result with a large deflection theory, 
Fig. 2. Although the two approaches are comparable 
in some ways, it was the buckling computation which 
had the more popular appeal. 

The picture began to change when Shanley® * intro- 
duced the concept of considering the loading process 
itself by returning to the strut problem and following 
what happens as the load is increased and the elastic 
limit is exceeded. He demonstrated conclusively that 
the classical buckling approach, which gives the critical 
load Px is not appropriate for the plastic range in which 
both nonlinear and effects appear. 
Perhaps, however, too much attention was paid to the 
remarkable proof that an initially perfect column could 
start to bend at the tangent modulus load, Pr. The 
buckling aspects of Shanley’s problem seemed to over- 
shadow the maximum load computation, or more gen- 
erally, his load-deflection relation, possibly because the 


nonconservative 
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STABILITY QF 


spread between P, and Px is so small for a column. 
As already mentioned, initial imperfections are always 
present so that basically the question of practical im- 
portance is how the deflection grows with load and 
what is the maximum load, P,, which can be carried. 
Pearson‘ contributed to the overall picture by proving 
that Shanley’s load-deflection curve for an initially 
straight column is the limiting load-deflection curve 
for an imperfect column as the deviation from straight- 
ness approaches zero, Fig. 3. 

In the Shanley example, it is nonlinearity which is 
responsible for deviation from straightness at the 
tangent modulus load, nonconservatism which pro- 
duces Px,. Nonlinearity is so prominent in the elastic- 
plastic range because a relatively small amount of de- 
formation produces the change from elastic response 
with its high modulus of elasticity to the plastic re- 
sponse at low modulus. 

Onat and Drucker® added the more elaborate ex- 
ample of a cruciform column in the plastic range which 
fails by twisting, Fig. 4. This plate problem solved 
by small deflection theory demonstrated the need to 
take into account what might be termed small but 
finite deformations which occur as the loading proceeds 
if extremely small initial imperfections are present. 
Here the spread between Py, and Px is large and of real 
practical significance similar in a way to the shell 
problem of Fig. 2. 

Hoff* and Ziegler” * focused attention on dynamic 
loading or dynamic disturbance superposed on static 
loading and pointed up the shortcomings of classical 
theory. A kinetic criterion of instability can be em- 
ployed to obtain correct answers for systems with per- 
fect geometry. The necessity or convenience of a 
kinetic criterion for essentially static problems, how- 
ever, remains to be considered. 


NON-LINEARITY, STATIC ANALYSIS 


As stated, plastic action is always both nonlinear and 
nonconservative. The first model to be considered, 
Fig. 5, has both characteristics but it is nonlinearity 
which provides the interesting features of its be- 
havior. Rod OA is rigid, the pin at O is frictionless 
and the spring attached to A has an elastic and a plastic 
range of force F shown schematically by the full line 
in Fig. 6. Small displacements only will be analyzed 
so that L@ replaces L sin 6, L cos @ is taken as L for the 
lever arm of F, and the position of A below its maxi- 
mum possible height is approximated by L6*/2. Upon 
unloading, the force-displacement curve is straight and 
parallel to the original elastic lineasshown. The actual 
force-displacement curve will be replaced by the broken 
line segments BC, CD, DD’, etc., for the convenience 
of description and of algebraic manipulation. Hoff® 
has previously made use of the model of Fig. 5 to dis- 
cuss ideal plasticity and creep. 

The mathematical expression of the simplified force- 
displacement relation depends upon whether the mate- 


INELASTIC SYSTEMS 


Si] 
wes 
si] 


rial is behaving elastically or plastically and is 


dF k,Lde (1) 


Il 


or 


dF kLdé (2) 


Expression (1) applies for FdF > 0 and 6 > @, on first 
> F, and 
Expression 


loading or more generally FdF > 0 and | F 
greater than any value reached previously. 
(2) applies for FdF < 0 or for |F| < F, or any pre- 
viously attained value. 

Suppose that the spring is so adjusted that F = 0 
when 6 = 0. The bar OA will then be in equilibrium 
in the vertical position for all values of load, P. The 
question of the stability of the equilibrium is not as 
trivial as might at first appear. If the usual static 
criterion is employed, OA is rotated through a very small 
angle @ and the work done by P, PL6*/2 is compared 
with the energy stored by the spring, FL@/2 where F = 
RLO. 


Equating the two terms leads to the critical load 
Pr = kL (3) 


The Shanley concept of increasing P as @ is introduced 
does not change the result at all. It is the stability of 
deflected positions which is the key to the physically 
significant behavior of the system. Equilibrium re- 


quires PL@ = FL or 
F = P9 (4) 
For the elastic range, 6 < 6,, F = kL@ and P = kL. 
When @ reaches or exceeds 6,, F = F, + k,L(@ — 4), 
where F, = kLA,. 
P§ = RL, + k,L(@ _ 6;) = (k — k,)L0, — 
k,L@ = kLO — (k — k,)L(@ — 6) (5 


as depicted in Fig. 7 by lines @ = 0. 
If the spring is so adjusted that F = Oat 6 = &, 


F = kL(@ — 0) = PO 


and 
A i. 8 
1— (P/kL) 1 — (P/Px) 


(6 


Beyond F,, where as shown dashed 
[A:/ (00 + )], 


in the elastic range. 
in Fig. 7a, P/Px = 


F = kLO, + k,L(6 — 6% — 6) = PO 


so that to maintain equilibrium P would have to de- 
crease as 9 increases. The maximum value of P there- 


fore is given by 
Py Px = (0; (A - 6;) | (7 


with the restriction Py, > k,L. 

In a sense Figs. 7a and 7b which differ only in hori- 
zontal scale are a graphical representation of the two 
extreme possibilities. If @, is large compared with 
probable initial and subsequent deflections P = kL is 


a valid buckling load and will be found experimentally. 
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On the contrary, if 6, is very small compared with 











hov 

probable initial deflection or equivalent eccentricity or P eri 
inclination of load P, the experimentally determined D 

critical load or the one found in practice will be close to 

k.L. Whenk,<k, the difference between the extremes | wh 

is very large. To this extent, the physical phenomenon 8 1] 

may be said to be strongly dependent upon initial im- L and 

perfection. In particular cases as for the cruciform, 2 on 

Fig. 4, it may turn out that the maximum load is _— 

bounded quite closely in practice. as . mar bar 

If the model is taken as important in itself instead of 8. 

simply a schematic illustration it is interesting to see G y ! 
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STABILITY OF 


how the extreme cases arise. F, = kL6#, may be re- 


written as 
0, = (F,/kRL) = (6./L) (8) 


where 6, is the maximum elastic elongation of the spring. 
A long soft spring will provide a 6, of several inches 
and Px will be obtained. Short snubbing bars or wires, 
on the other hand, permit a 6, of the order of thou- 
sandths of an inch and unavoidable measures in the 
bar and the loading will give an effective 6) many times 
§,. Py will then tend to be close to &,L. 

Although Fig. 6 shows the unloading behavior to be 
nonconservative, nothing in this section involved un- 
loading so that a nonlinear elastic material would 
exhibit exactly the same effects. 


NONLINEARITY, KINETIC ANALYSIS 


The model of Fig. 5 will now be considered from the 
alternative kinetic point of view. Now the question 
to be asked is whether a dynamic disturbance applied 
to a perfectly aligned system in equilibrium produces 
bounded or unbounded oscillations within the frame- 
work of small deflection theory. The previous static 
analysis makes it quite clear that the answer will again 
depend upon the magnitude of the disturbance just as 
the size of 6) was significant. 

Suppose that under constant load P, with @ = 0 
the bar OA is given an initial velocity #. Conserva- 
tion of energy gives, for small 6, 


1(6? — 6°) = PL6?/2 — kL@?/2 (9) 


for? < @, 
and 
I(62 — 6?) = PLe@?/2 — 


RLO,7/2 + RLO, (6 — 0:) + k,L(@ — 6,)?/2] (10) 





for@ > @,and6> 0. Jis the moment of inertia of the 
bar OA about 0. 

The system is unstable, that is, motion in the initial 
direction will not stop, 6 > 0, if the energy stored and 


dissipated in the spring cannot be as large as 
16° + PL@/2 (11) 


Obviously, from Eqs. (9) and (10) P < Px = kLisa 
requirement for stability. At P = Px, the slightest 
% will cause collapse, at P > Px the system runs away 
quickly. A plot of @ vs. 6 is most convenient for ex- 
hibiting such properties, Fig. 8. 

If P < k,L, there is no 6)? which cannot be absorbed 
in accordance with small deflection Eq. (10). The 
system is therefore stable for P < k,L and the 86, 6 dia- 
gram is composed of portions of ellipses, Fig. 9. For 
kL <P < kL there will be levels of initial disturbance 
above which the system will not recover. Below these 
values of 6) the bar will come to rest and return. It 
is clear on physical grounds that on the return motion 
the system will again stop and oscillate back and forth. 
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Unless the system collapses on the first try it is stable, 
Fig. 10. 
gous cases when the external forces applied are con- 


It is for this reason that in this and analo- 


servative, static and kinetic analyses are equivalent. 
Fig. 10 also shows the effect of combinations of initial 
velocity and initial displacement. 


NONCONSERVATIVE ASPECTS 


Fig. 11 shows the model employed by Duberg and 
Wilder,® following closely the model and concept of 
Shanley, to investigate the effect of initial imperfec 
tion in columns. DLGR is rigid and G is constrained 
to move vertically only by means of a frictionless guide 
so that the system has two degrees of freedom y, @. 
L and R rest on two identical short bars of work- 
hardening material whose force-displacement dia 
grams are each idealized as two straight lines, Fig. 12. 
If a bar has been loaded to B and loading continues 


F = Fg + k,(6 — 6,) (12 
while if load 1s removed 
F = F, — k(ép — 5) (13) 


The static analysis is well known and basically has 
already been described in Fig. 3 where 6 should be re- 
placed by 6. The load rotation relation for the simpli- 
fied model is given in Fig. 13. An initially imperfect 
system, 6 = 4, will rotate more as load is applied and 
will have deviated appreciably from the vertical by the 
time the tangent modulus load, P7 (Fig. 14) is reached; 


Pr = k,b?/L (14 


A perfect system can deviate from the upright position 
at P, and the slope of the P, 6 line will be progressively 
flatter the larger the value of P at which rotation is per- 
mitted (see line segments between P; and Px of Fig. 3). 
At the reduced modulus load P = Px, Fig. 14, the slope 
would be zero; 


Px = (k,b?/L) [2k/(k: + B)] (15 


The third critical load of interest is the Euler elastic 
buckling load which can be obtained by putting k, = 
kin Eqs. (14) or (15); 


Py = kb? L (16) 


The preceding discussion has implicitly assumed that 
P, and therefore Px and Py are above the yield force 
2F, (Fig. 12). Should P, be less than 2F,, the initial 
imperfections may be small enough so that the mate- 
rial may not be made aware of the plastic deformation 
it could experience. Nothing unusual will then happen 
at the tangent modulus load. 
is a perfectly plastic or non-work-hardening material, 
k, = 0, for which P; = 0. In what follows, therefore 
Pr > 2F,,. 


A dynamic investigation of the geometrically perfect 


The obvious example 


model is more interesting than the well-explored static 


study. As in the analysis of the first model, the con- 
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sideration of a kinetic criterion of stability leads to in- 
vestigation of the motion which follows arbitrary initial 
disturbances, and the question to be answered is 
whether or not the ensuing motion is limited. The 
system has two degrees of freedom, one of motion of 
G vertically, y, and one of rotation, 0, and the arbi- 
trary initial dynamic disturbance will be characterized 
by wo, 6 where dot indicates time derivative. This 
pair of values determines the initial instantaneous 
center of rotation C and of course the initial angular ve- 
locity of the rigid bar DLGR. If the x-coordinate of 
C is between L and R, one supporting bar will start to 
unload and the other to load. On the other hand, if 
the x-coordinate of C lies outside L to R both bars will 
start to load or unload depending upon the sign of 6. 

It is convenient in the analysis to take the center of 
mass at G. The equations of motion for small dis- 
placements y, 0 from the untilted configuration under 


load P are then 


my = —(kt + kre)y + (Re — kx)b0/2 


(17) 


16 = PLO/2 — (ki + kp)b?0/4 + (kr — kx)yb/2 (18) 
where the bar constants are denoted by k, and kp, each 
of which may be either k or k,. The model therefore 
behaves like a linear elastic system at the beginning of 
the motion except that the nonconservative stress- 
strain relations determine the proper value of the k’s. 
The kinematic requirement is that k; = kr unless 
y| < [b0/2). 

From Eq. (17), positive vertical acceleration j re- 
quires 


: (Rr — k,)/(ki + kr) | 60/2 (19) 
and from Eq. (18) positive angular acceleration 6 re- 
quires 


P > (kr + kr)b?/2L — (kr — kr)by/Le (20) 
where it is assumed that 0 > 0. 

The smallest value of P for which the angular ve- 
locity will not decrease as the motion begins is therefore 
seen to be the tangent modulus load P;. This value is 
obtained from Eq. (20) by taking both LZ and R to 
shorten so that k; = kre = k; However, under these in- 
itial conditions the vertical acceleration will be negative 
and the motion will reverse because bar R must even- 
tually unload. A contained oscillation will thus be 
set up if P is less than the reduced modulus load Px. 
Above Px, positive angular acceleration is accompanied 
by positive vertical acceleration as one bar loads, k, = 
k,, and the other unloads, kr = k. 

It is interesting to note that the elastic Euler load, 
Px, may also be obtained from Eq. (20) by assuming 
an initial disturbance that extends both bars L and R; 
ki = ke = k. Positive angular acceleration for these 
initial conditions requires P > Py. Of course, the 
model will collapse for any P > Px because the upward 
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Fic. 14. Tangent and reduced modulus load computation, 


motion will reverse and failure occur in the conven- 
tional manner. 

A discussion of the behavior of the model in terms of 
energy is also useful. Designating the change in kinetic 
energy for the very beginning of the motion by A7, 
first order terms cancel out and 


AT = PL6?/4 — kily + 60/2)?/2 — 
Rr (y — 00/2)?/2 (21 


The total kinetic energy will therefore decrease ini 
tially for P < Px no matter what the initial disturb- 
ance, although the rotational kinetic energy increases 
initially for y > 66/2 and P > Py. Above Px, both 
the rotational and translational energy increase for 
y < 00/2. In summary, therefore, the motion of the 
system must be examined under all possible disturbances 
to understand the significance of the critical loads ob- 
tained with a kinetic criterion of stability. For the 
model of Fig. 11, because k, is constant the motion 
settles down when Py < P < Px whether the incre- 
ment in rotational energy is initially positive or nega- 
tive. With the usual curved stress-strain diagram it is 
clear that collapse would occur for P less than Px by an 
amount depending upon the magnitude of the initial 
values 6, ¥. Nevertheless it is interesting and im- 
portant to note that all three critical loads, P;, Px, and 
P, can be obtained from an analysis of a geometrically 
perfect subjected to infinitesimal dynamic 
disturbance. 


system 
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The Estimation of Normal-Force, Drag, 
and Pitching-Moment Coefficients for 
Blunt-Based Bodies of Revolution at 

Large Angles of Attack 


HOWARD R. KELLY* 
U.S. Naval Ordnance Test Station, Inyokern 


ABSTRACT 


A method is developed which accurately predicts for blunt 
based bodies of revolution the normal-force coefficient and 
pitching-moment coefficient for angles of attack far beyond the 


range of potential theory. The method is based on the super- 


position of the results of potential theory and the viscous force 
on a cylindrical body due to the transverse component of flow 
In contrast to previously used methods, the viscous cross-force 
is assumed not to be in a steady state but in a transient develop- 


ment along the body 

This method is compared with experimental data for both 
subsonic and supersonic flows and for both laminar and turbulent 
boundary layers. An extension of the method to the prediction 
of the drag-coefficient increment due to yaw is made and is 
compared with subsonic data. The deviations between theory 
and experiment are in most cases an order of magnitude smaller 


than for previously used methods 


INTRODUCTION 


1. AERODYNAMIC DESIGN of aircraft rockets has 
some unique features in that the most practical 
body shape from other design considerations usually 
isa cylindrical shape with a blunt base and a moderately 
blunt nose. Little good theoretical information is 
available on the lift and pitching moment for bodies 
of this type, except at extremely small pitching angles. 
For this reason, a considerable number of wind-tunnel 
tests have been made by the Naval Ordnance Test 
Station to build up a fund of information for design 
purposes. 

While examining some of these data for the purpose 
of correlation between tests and with available theory, 
it was found that certain subsonic test data on short 
bodies at large angles of attack differed widely from 
any theoretical predictions available. The theory of 
Max Munk! was found to give good results at small 
angles; this is satisfactory, since potential theory is 
not expected to be valid for large angles. A theory 
proposed by H. Julian Allen as early as 1949, and 
published in reference 2, which has enjoyed a measure 
of success for supersonic flow data, grossly overestimates 
the effects of viscous cross flow at low subsonic speeds, 
where it would at first be expected to be at its best. 


Presented at the Aerodynamics Session, Twenty-Second 
Annual Meeting, IAS, New York, January 25-29, 1954 


Aerodynamics Branch 


The theory of reference 2 has been carefully ex 
amined, and a completely new version of this theory 
is presented here. The problem of viscous cross flow 
has been reformulated with quite different results. 
The use of Munk’s results for the linear potential theory 
contribution is a fairly good approximation for long, 
For the bodies considered here, 
The results 


slender head shapes. 
a revision of this policy is recommended. 
of wind-tunnel tests by NOTS are presented in support 
of this new theory, showing excellent correlation. 
The predictions of Allen’s theory are included for 
comparison, and the results of Munk’s linear theory 
are included as a sort of reference line, without intend 
ing a comparison. 

This method is found to give formulas for certain 
cases which would appear to make it more difficult 
to use for engineering design purposes than Allen's 
method. It is shown, however, that considerable 
simplification can often be effected for this type of 
application, without seriously impairing the accuracy. 
An important result of this method is that it provides 
a means of extrapolating the results of tests at low 
Reynolds Numbers, with completely laminar boundary 
layer, to high Reynolds Numbers, with fully developed 
turbulent boundary layer. This is a valuable asset 
for design work when test facilities are limited. 


HISTORICAL DEVELOPMENT 


Until recently, there had been little theoretical 
work done on theories of lift for bodies of revolution, 
except the linearized potential theories. A notable 
early example of potential theory was that of Munk,’ 
which is useful for incompressible flows at small angles 


of attack and is often approximated by the expression 
Cr, = 2a (] 


for slender bodies without boat-tail. 

Munk realized that his linearized potential theory 
would apply only at small angles of attack and that 
it would fail at large angles, chiefly because of viscous 
There was apparently little interest at that 


effects. 
Recently, 


time in developing a theory for large angles. 
Allen? has developed a method for predicting lift, 
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pitching moment, and the drag increment due to yaw 
which has enjoyed considerable success with slender 
bodies of revolution at supersonic speeds. 

Allen has assumed that the viscous contribution 
to lift is independent of the contribution of the potential 
flow and writes the lift coefficient and pitching-moment 
coefficient as the sum of two terms, 


SB a A 
’ “ ° . 9 s+ «@ 
Cr = (ko — ky) A sin 2a@ cos = + nCb- 4 sin? @ COS @ 


(2) 


V — Sp(l — 
Cu = | = 


Xm) q a 
F - | (ko — k,) sin 2a@ cos : oF 


. A, 
nC De A 


(Xm — X,) sin? a@ (3) 


where the first term is always the result of Munk’s 
theory, modified by the results of Ward,* who showed 
that the potential cross force is directed midway 
between the normal to the axis of revolution and the 
normal to the wind direction The quantity Cp, as 
used by Allen refers to the steady-state drag coefficient 
of an infinite cylinder in a flow normal to its axis. 
This was found by him to be Cp, = 1.2 for a wide range 
of Reynolds Numbers corresponding to the usual 
range of cross-flow Reynolds Numbers being con- 
sidered. The factor n corrects this cylinder cross-flow 
drag for the finite fineness ratio. Values of 7 have 
‘been tabulated by Goldstein.‘ The x,, is the (arbi- 
trary) point about which pitching moments are meas- 
ured; x, is the axial position of the centroid of the 
plan-form area A,; and the factor sin* a arises from 
the transformation of dynamic pressures from the 
cross velocity U sin a to the free-stream velocity V. 
Since his theory is only approximate, Allen simplified 
his formulas by making the approximations sin a = a, 


ko — ki = 1, and cos a = 1. This results in 


S 4 
Cr = 2 ( *) a+ Cp, (“”) a? (4) 
y 
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Geometry for body of revolution with angle of attack a. 
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V — Spall — x,,) 
Cu = 2/ “ : ‘Jat 


‘ A, 
nCp- { (Xm — X-)a~ (5 


and his theory is used in essentially this form. J. A. 
F. Hill, at the Naval Supersonic Laboratory, Massa- 
chusetts Institute of Technology, used Allen's results 
in 1950° but modified them to include boundary-layer 
displacement thickness in the potential theory term, 
thus changing the effective values of both Sg and J’ 
He also observed that the assumed value of Cp, = 1.2 
is not appropriate if the axial boundary layer is turbu- 
lent, even though the cross-flow Reynolds Number 
corresponds to a laminar boundary layer. 


THE Viscous Cross FLow 


In the theory presented in this paper we shall be 
chiefly concerned with the normal-force coefficient 
Cy rather than the lift coefficient C,. This will 
eliminate the factor of cos a in the second (viscous 
term of Eq. (2). It also eliminates the necessity of 
Allen’s assumption that the component of axial drag 
involved in the lift force may be neglected. This 
assumption would be invalid for the angular range of 
interest here. As Hill has pointed out, Allen’s Eq. 
(4) would be more appropriate for the normal-force 
coefficient Cy. 

Another modification of Eq. (4) is evident from a 
Munk’s theory. Munk 
‘as such but derives 


closer examination of does 
not once mention the “‘lift force’ 
the expression for the pitching moment, which is a 
couple for a closed body. Although the resultant 
force may be midway between the normals to the axis 
and wind direction, as shown by Ward, any axial com- 
ponent will not be involved in an expression for normal 
force or in Munk’s expression for pitching moment. 
Therefore, the factor cos a/2 is not appropriate for 
the normal-force coefficient. A simple resolution of 
forces shows that it is also inappropriate for the hit- 
force coefficient. 

Allen’s fundamental equation for the viscous cross 
flow is useful if properly interpreted. He writes the 
viscous contribution to the cross force from a cylindrical 
element of length dx as 


dF = 2rCp.(pV,7/2) dx (6 


where 7 is the body radius at the point x and J’, is the 


cross velocity, as shown in Fig. 1. The quantity p is 


the usual mass density, and Cp, is the drag coefficient 
of a circular cylinder at the Reynolds Number 


Re, = 2rV./9 (7a 
and the Mach Number 
M. = V./a (7b) 


with the kinematic viscosity v and speed of sound a 
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Approximation of Schwabe’s results by a fifth-degree 
polynomial 


evaluated in the undisturbed flow. Then, since 


Ve. = U 


sin a (8) 
it follows that 
dF = 2rCpq sin® a dx (9) 


in terms of the dynamic pressure g = pl*/2 of the free 
stream. 

Since the drag coefficient of a circular cylinder in a 
cross flow is usually given for the two-dimensional case 
of an infinite cylinder, Allen adds a parameter 7, 
which represents the reduction in drag coefficient of a 
cylinder due to its finite length. Then, if the moment 
coefficient for a body of diameter d, referred to a point 


Xm 


Cy — M, gAd (10) 
and the normal-force coefficient 
Cy = N/qA (11) 


we find the increments in these coefficients, due to 


VIsc( sity, to be 


: 2nsin?a {! _ 
ACy = 1 rCp- dx (12) 
wet 7] 
: <7 Sin” @ " - 
ACy = rCpo-(Xm — x) dx (13) 
Ad 0 


which lead to the final terms of Eqs. (4) and (5). 
These terms are entirely due to viscous effects, since 


Cp, always vanishes for a true potential flow. 


EVALUATION OF PARAMETERS 


The departure of the theory of this report from that 
of Allen lies principally in the evaluation of the pa- 
rameters of Eq. (12) and of the potential theory. It was 
well known to both Allen and Hill that the coefficient 
Cp. should not be the steady-state drag coefficient, 
as used by them, but should be related to the transient 
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effect found by Schwabe.® The results of Schwabe's 
experiments are reproduced here in Fig. 2. The solid 
circular 


curve represents the drag coefficient for a 
cylinder moving crosswise in a fluid when started 


impulsively from rest, as computed from the observed 
flows. 
drag coefficient rises rapidly from zero until the pa 


For a cylinder of radius 7;, it is found that the 


rameter I|,f/r; is about 4. <A slow increase continues 
until when this parameter reaches a value of 9, the 
limit of validity as shown in Fig. 3, the drag coefficient 
is 2.07 or twice the steady-state value of about 1.0 for 
the Reynolds Number of the experiment. 

Referring to Fig. 1, we may consider a plane lamina 
of air of thickness dx moving along the body axis 
with speed l’ cos @ and across the body with speed 
Usina. At the nose of the body, this lamina suddenly 
begins to flow about a circular cylinder of variable 
A potential cross flow will exist at the nose, 
flow 


radius r. 
with a gradual development of 
Schwabe’s parameter becomes 


viscous cross 


along the body. 
V. U' sin a x x 


t =— - — 
r r U cos a r 


tan a@ (14 


In order that Eqs. (12) and (13) may be readily 
integrated, it is advisable to approximate Schwabe’s 
results by a simple function. The experimental curve 
by Schwabe resembles a portion of an are tangent 
function; 


curve, which suggests the use of such a 


however, the evaluation of the integrals obtained 
appears to yield more complicated formulas than the 
present problem deserves. The Taylor expansion of 
an are tangent function 


y = tan! x = x — (x*/3) + (x5/5) + (15) 


suggests the use of an odd polynomial, and the factors 
containing tan @ separate out immediately this way. 
The limits of applicability of this theory, due to a 
change in flow conditions at large angles, renders the 
use of terms beyond the fifth power questionable, so 
that a polynomial of the form 





FINENESS RATIO, £ 
ae T 





20} --——++ 

















| 
| 
0 i 
0 5 10 15 20 25 30 35 40 45 
ANGLE OF ATTACK, @, DEGREES 


Fic. 3. Limit of validity of present theory 








552 JOURNAL OF THE AERONAUTICAL 
; _s ; x3 : 
Cpe = A tana+ B’ |. tan'a + 
, si 


Cc’ tanba + 


5 


(16) 


will be sufficient for practical purposes. 

The cylinder cross-flow drag function may be inte- 
grated for specific body shapes. Since the integral is 
multiplied by sin? a, we may simplify the results bv 


using 
sin? atana = a’®+ (a7/15)+... (17a) 
sin? a tan? a = a’ + (2a7/3) + (17b) 
sin’ atanha = ai + (17c) 


so that, to terms in a’, the normal-force and moment 
coefficient increments will be of the form 


Atty = (2Cp,- A [Dy ai + Kya’ + Fy a‘ ] (18) 


ACy = (2Cp,/A)(Dua’ + Eva’ + Fya’| (19) 
where Cp, is the steady-state value of the cylinder 
cross-flow drag coefficient. 

It is also necessary to choose an appropriate value of 
the parameter yn for the transient phenomenon being 
considered. The fact that 7 is different from unity 
for the steady state is principally due to the distortion 
of flow near the ends of a cylindrical body. A lamina 
of fluid will not remain in a plane but will flow in such a 
way that the total flow will meet less resistance. In 
the case of the transient development being considered 
for a blunt-based body of revolution, this effect will 
not have time to develop at the base before the flow 
leaves the body. A distortion may develop at the 
nose, but the viscous force contribution at the nose is 
extremely small, and the effect may usually be neg- 
lected. For these reasons, the factor 9 will be replaced 
by unity. 

If a body does not have a blunt base 
boat-tail—the above argument will fail, since con- 
siderable distortion of the cross-flow planes will occur 
in a region where the viscous effects are greatest. For 
this case, the value of » may differ significantly from 
unity, which is one reason the present report is re- 
stricted to blunt-based bodies. 

There remains the problem of evaluating or replacing 
(ko — k,), the factor in Munk’s theory representing 
the difference in apparent mass for longitudinal and 
transverse motion of the body. It has been suggested 
by Prof. Ralph Upson, and confirmed by Bernard 
Leadon, both of the University of Minnesota, that an 
analysis similar to that of Upson and Klikoff’ will be 
more appropriate than Munk’s for half-bodies of the 
type considered here. In brief, the factor (k2 — kj) 
should be replaced by (AB/2) cos? 6, where A = 
1+ kh, B= 1+ k, and tan 6 = dr/dx. The correction 
factor will thus be variable along the nose of the body, 
except for conical nose shapes. For extremely long, 


that is, has 
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slender nose shapes, this new factor may be replaced 
by unity without serious error, just as could be done 
with Munk’s (k» — k,). 


BOUNDARY-LAYER CORRECTIONS 


Hill has pointed out that the boundary layer may 
necessitate two corrections to the theory of Allen 
one in the viscous flow term and one in the potential 
term. His modified theory does not eliminate the 
fundamental difficulty and principal source of error, 
but this type of correction may not always be neglected. 

Allen’s viscous cross-flow term assumes that the cross 
flow is completely independent of the axial flow and 
will therefore depend only on the cross-flow Reynolds 
Number [Eq. (7a)] and Mach Number [Eq. (7b 
There is good experimental evidence, however (ref- 
erence 4, p. 431), that the critical Reynolds Number 
for the decrease of the drag coefficient of a cvlinder 
from the laminar Jow value to the turbulent flow value 
changes rapidly with the turbulence level present in 
the streain. It seems logical that, if a cross flow in 
this range of Reynolds Numbers encounters a turbulent 
axial boundary layer, the drag coefficient 
Cp, to use will be that for a turbulent cross flow. The 
coefficient appears to change by a factor of about 3, 
and the steady-state values of Cp, that will be used 
here will be 1.2 for the case of completely laminar 
boundary layer and 0.35 for a (almost) completely 


proper 


turbulent boundary layer. 

The effect of boundary layer on the potential theory 
term arises from the fact that the flow is displaced as 
though the body shape were changed. In particular, 
the effective base area Sz and volume IJ’ will be greater. 
The important boundary-layer parameter will be the 
displacement thickness 

6* = fi’ [1 — (u/U)] dy (20) 
For laminar boundary layer this thickness may be 
found approximately from 
/ 
5*/x = A/~V/Re 2] 
where 
A = 1.73 + 0.48M? 22 


The turbulent boundary-layer displacement thickness, 
as shown by Prandtl from the seventh power law 


u/U = (y/6)'” (23 
is given by 

8*/x = K/(Re)°” (24 
where AK = 0.046 for incompressible flow on a flat 
plate. Flat-plate data from the Daingerfield super- 


sonic wind tunnel indicate an approximate value of 
K = 0.08 for a Mach Number of 1.56, which is of 
interest in this report. Hill has computed the frac- 
tional increase in base area from Sz to Sz’ and ap- 
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PITCHING-MOMENT 


proximate formulas for the fractional increase in 
volume from V to 1’, based on the assumption of 
constant diameter and fully turbulent boundary layer. 


These are, for the laminar case, 


S A / 


A l 
—l=4 1+ (25 
Sp AY Re, V7 | \ Re, A 


y’ 2nd* A si 3 A / 
—l= : ( L \ (26) 
Y’ 3 V V Rea \d 4~—V/ Re, Vd 


for the turbulent case 


Sp’ a x (‘)" [ | K 
Sr 7 ( Rea ).2 d bal ( Re. ho 


y’ ; 2 =) K (3) x 

yt~” | V / (Req)®:? \d 
9 K 1\o§ 

[ + : ( ) | (28) 
13 (Rez): \d 


SUPERSONIC POTENTIAL FLOW 


The proposition that the normal force and pitching 
moment may be found by the simple addition of the 
potential flow solution and the viscous cross-flow 
solution is an important contribution from the theory 
of Allen. The logical application of this principle 
will use the best available potential flow solution. 
For subsonic flow, the theory of Munk, corrected for 
boundary layer, gives excellent results. For super- 
sonic flow, however, Munk’s theory is not exact, except 
possibly in the limiting case of bodies with extremely 
long, slender noses. 


A much better approximation to the potential flow 
solution for supersonic flows is obtained by the second- 
order theory of Van Dyke.* This theory has been 
applied recently to the calculation of pressure dis- 
coefficients, and _ pitching- 


tributions, normal-force 


moment coefficients for a series of tangent-ogive- 
cvlinder combinations.’ 
computed for a finite number of discrete supersonic 
Mach Numbers, so a method of interpolation is de- 
sirable. In Fig. 4 the values of dCy/da for a = 0 


degrees as computed in reference 9 are plotted as a 


These must necessarily be 


V M? — 1, with ogival head length 
These values were computed for a 


function of B = 
as a parameter. 
total body length of 14 calibers, and do not change 
appreciably with increase in body length beyond 9 
calibers. No particular significance has been attached 
to the linear relation, except for a convenient method 


of interpolation. 


NORMAL FORCE AND PITCHING MOMENT 


The foregoing discussion may be summed up by 


writing the complete equations for normal-force 


coefficient and coefficient of pitching moment as 


COEFFICIENTS 


FOR BLUNT-BASED BODIES 553 





NORMAL-FORCE -COEFFICIENT SLOPE, Cy, 
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Fic. 4. Initial slope of coefficient of normal force from Van Dyke 
Theory, for / = 14 

: dCy j 

Cy = sin @ COS @ +> 

da aut 
yi 
- De . P . . 
i [Dya + kya’? + Fya’ (29 


sin a cos a + 


— (<=) 
> da 


a=( 
2¢ D 


A 


|[D ya? + Eya' + Fyai (30) 
The dCy/da and dCy/da are given by Munk’s theory 
at subsonic speeds, after correcting for the apparent 
mass factor and the boundary-layer thickness; for 
supersonic speeds they are found from Van Dyke's 
second-order theory, with boundary-layer corrections 
added. The D, E, F coefficients are computed from 
the integrals in Eqs. (12) and (13), and Cp, is 1.2 or 
0.35 as the axial boundary layer is laminar or turbulent, 
respectively. 

As an example of the results of the integration, the 
integral for the normal-force coefficient for a cone- 
cylinder is 


. 


[? 
/ (Code = A’ 5 tana + B'tl4 + h*] tan? a + 
te ; 
16C€ 6 +. = tan°’a (31) 
) oO 


where / and h are the fineness ratios of the entire shape 
and of the head, respectively. The polynomial coeffi- 
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cients have been evaluated as 


A’ = 0.49 
B’ = —0.0056 (32) 
C" = 8) 00003) 


Other examples and more detail of the integrations 
will be found in reference 10 and will be given ina 
forthcoming NAVORD report that will supersede 
reference 10. 

The equations become particularly simple for the 
case of subsonic flow with a < 10° and / < 15 calibers. 


Then 
Cy = 2a + 0.49(Cp,l?/A) a? (33) 
Cuz = (2V/A)a + 0.49(Cp./?/3A ) a? (34) 


which equations are as simple in form as Eqs. (4) and 

5) and agree much better with experiment in the 
range of validity mentioned above. 

Comparisons of the complete seventh-degree poly- 
nomial with subsonic data are shown in Figs. 5 and 6. 
The effects of boundary-layer thickness and of the 
apparent mass factor very nearly cancel and may be 
neglected for this case. Reference lines have been 
drawn corresponding to the linear terms of Munk’s 
theory to emphasize the nonlinear effects. Further 
comparisons of normal-force coefficient are made for 
supersonic data, chosen to represent turbulent and 
laminar axial boundary layers (Figs. 7 and 8). For 
the turbulent layer, Allen's formula appears to agree 
with the data. This is somewhat fortuitous, since the 
error from using Munk’s formula for the potential 
term nearly balances the error in the viscous term. It 
may be noticed that his predicted initial slope dis- 
agrees with the data. 

An additional important use of the present method 
is in extrapolation of experimental data from small 
models to large-scale models at large angles of attack. 
The initial terms of Eqs. (29) and (30) may well be 
wind-tunnel data, corrected for 


obtained from 


boundary-layer conditions. 


ESTIMATION OF DRAG INCREMENT 


The extension of this method to the prediction of 
wind-axis forces follows immediately from the pre- 
ceding discussion but is less straightforward in that a 
knowledge of the drag at zero angle of attack is re- 
quired. Resolution of force vectors for the geometry 
of Fig. 1 shows that, for subsonic flow, 


— Cp,(1 — cos* a) + 


2 sin’ @ sas ™ 
1 rCp.dx (35) 


The factor tan a/2 is due to the result of Ward that 
the potential force vector is at an angle of a/2 from the 
The final (viscous) term, 


’ 3 ‘ a 
Alp = sin 2a tan 


normal to the wind stream. 
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when integrated, yields an expression of the form 
20s. 


1 [Dp at + Epa’ + Fopa’!] (36) 


(ACp)y = 
so that the result will be a polynomial in even powers 
of a. Considering the special case of subsonic flow 
with a < 10°, / < 15, the term in Cp, usually becomes 
small, and one may use the approximate result for 
blunt-based bodies, where Sz = A, 


A )l? a! (37) 


ACp = a? + 0.49(Cp, 


which is as simple to use as the corresponding equation 


from the method of reference 2, 


ACp = a? + nCp.(A, A)a (38) 


There does not appear to be a significant difference in 
the accuracy of those equations, but the distinction 
becomes apparent at values of 2/ tan a > 5, where the 
complete polynomial must be used, as shown in Fig. 9. 
This method has not yet been tried for supersonic 
drag increment, but the extension to this case is fairly 
obvious. 


CONCLUSION 


A method has been developed for the estimation of 
normal-force coefficients, drag coefficient increments, 
blunt-based 


and pitching-moment coefficients for 


bodies of revolution at large yaws. The equations 
reduce to 
Cy = 2a + 0.49(Cp,/?/A ) a’ 
Cuz = (2V/A)a + 0.49(Cp,/*/3A ) a? 


ACp => a’ + 0.49( Cp, A \/*q@? 


} 


for subsonic speeds if a < 10°, 2/ tan a 
as the boundary-layer thickness does not increase the 


5, as long 


base area appreciably. 

It is important to notice that a seventh power poly- 
nomial is recommended if 2/ tan a > 5, and that, for 
the case of supersonic speeds, a more accurate potential 
theory or experimental data should be used for the 
potential contribution. Boundary-layer thickness may 
be important in some cases. 

An important result of this method is that it provides 
a means of extrapolating test data from small models 
with laminar boundary layers to free-flight conditions 
with turbulent boundary layers. The method is 
limited to bodies without boat-tail, but this happens 
to be a practical shape in the design of many types of 


aircrait rockets. 
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The Static Stability of a Two-Dimensional 
Curved Panel in a Supersonic Flow, 
with an Application to Panel Flutter 
Y. C. FUNG* 


California Institute of Technology 


ABSTRACT oe = dimensionless coefficient for initial deflection, see Eq 
(12). A; indicates amplitude of sinusoidal arch 


An elementary analysis is given for the problem of steady-state . = ratio (change of axial thrust)/(Euler load), see Eq 
stability, in a supersonic flow, of a two-dimensional curved panel, (28). 
hinged at two opposite edges, with the distance between the : 
hinges held fixed. The variation of the panel shape with in- > = indicates summation over m with the term » = m 
creasing dynamic pressure of the flow is studied. The existence of n= 1 deleted 


a critical supersonic speed of flow, above which it is impossible for 
a bent panel to maintain static equilibrium, is confirmed (1) INTRODUCTION 
Application of the solution to the panel flutter of a flat panel of 

supersonic aircraft due to thermal buckling is discussed. The vies PROBLEM TO BE studied, in its simplest form, is 
critical dynamic pressure of static instability is an upper bound to as follows: Consider a long strip of thin plate of 
that of panel flutter. Since this upper bound approaches zero at span L,a portion of which is shown as ABCD in Fig. 1. 
the initial stages of thermal buckling, it is concluded that if 
panel flutter is to be avoided, then, in general, thermal buckling 
of a flat plate cannot be tolerated. 


The plate is made of uniform, perfectly elastic material 
and has uniform thickness. Assume that the condition 
is two-dimensional so that the plate can be represented 
by a section of unit width such as EF in Fig. 1. The 















: F i ord haven a decks ee -_ leading and trailing edges of the plate (edges AB and 
am, bm = Fourier coefficients describing the initial and final : ; Edie 
plate deflections, respectively CD respectively) are hinged to rigid supports so that 
ee = Fourier coeflicients describing external loading the span Z remains constant at all times. Initially the 
p = external load per unit area 
Pa = aerodynamic pressure 
q = dynamic pressure = (1/2)pl? B 
t = thickness of plate 
¥, Yo = actual and initial curve of the mid-surface, respec- U we C 
tively | 
B. = dimensionless coefficient for final deflection, see Eq. 
(12) 4 
D = Et#/12(1 — yu?) F 
E = Young's modulus 
H = axial thrust per unit length 
Ho = initial axial thrust per unit length 
L = width of plate, span of arch 
M = Mach Number (a) 
Q = dimensionless coefficient for aerodynamic loading, 
see Eq. (12) 
Qe = critical value of Q above which the plate becomes 
unstable 
Qer = an upper bound of Qe y 
R = dimensionless coefficient for external load, see Eq. 
(12) 
5 = dimensionless coefficient for axial thrust, see Eq. 
(12). S = 1 implies the fundamental Euler 
buckling load 
B = dimensionless coefficient of support rigidity with re . 
spect to axial displacement, see Eq. (10). 8 = 1 
implies infinite rigidity 
A = axial displacement of supports & _— — 
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STABILITY OF CURVED 


plate is flat and free from stress. It is then uniformly 
heated so that an axial thrust (compressive) of intensity 
H{, per unit length is introduced into the plate. If //, is 
less than the Euler buckling load of the plate, the plate 
If [7p is equal to the Euler buckling load, 
The heated 


remains flat. 
the plate buckles into a sinusoidal arch. 
plate is then exposed to a steady supersonic flow parallel 
to the unbuckled plate. The plate deflects under t! 
aerodynamic pressure. The problem is to determine 
the stable deflected configurations of the plate. 

Since only static equilibrium is considered, the order 
of events has no effect on the final outcome. Hence 
the problem named above is the same as the case in 
which the heating of the plate is caused by the flow 
itself. Furthermore the thermal strain may be replaced 
by mechanical strain such as caused by a displacement 
of one of the supports. 

The center of interest lies in the stability of the plate. 
The possibility of instability arises from the nature of 
supersonic aerodynamics. Let us assume that the de- 
flection is so small that the aerodynamic force can be ob- 
with sufficient from the linearized 
In a two-dimensional supersonic flow the aero- 


tained accuracy 
theory. 
dynamic pressure is directly proportional to the slope of 
the solid surface with respect to the undisturbed flow. 
For a symmetric arch, the supersonic aerodynamic load- 
ing is thus antisymmetric, tending to deform the plate 
into a more complicated configuration. It is not obvious 
that such an aerodynamic pressure distribution is de- 
stabilizing, but it will be shown that when the dynamic 
pressure g increases, there exists a limit beyond which a 
buckled plate with fixed supports simply cannot main- 
tain any equilibrium configuration. 

The problem so formulated can be applied to an aero- 
elastic instability called “panel flutter,” which was ob- 
served in earlier rocket tests and has been confirmed re- 
cently in controlled wind tunnel experiments. The thin 
skin of a wing or fuselage structure, when exposed to a 
supersonic stream, may develop various types of oscilla- 
tions, which may cause fatigue and destruction. 

One type of panel flutter arises in practice from the 
bending or buckling of the skin due to aerodynamic 
heating. The behavior of a bent skin panel in a super- 
sonic stream depends critically on the nature of the sup- 
porting framework. If the edges of the panel are fixed in 
space, the membrane stress (i.e., the thrust in the skin) 
varies with the amplitude of the deflection and the 
problem is nonlinear. If the supporting structure is 
‘perfectly flexible,’ so that no resistance is offered to 
displacements of the panel edges in the plane of the 
skin, the membrane stress is independent of the ampli- 
tude of the deflection and the system is linear with re- 
spect to aerodynamic loading. 

As an introduction to the study of panel flutter of the 
“rigid-support”’ case, the static equilibrium of a two- 
dimensional curved panel may be considered. The 
argument is that when the dynamic pressure q is so high 
that no stable static equilibrium is possible, motion en- 
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sues. This motion (panel flutter) cannot be a stable 
one; for if it were, then after the elapse of a finite time 
interval a static equilibrium configuration will prevail, 
which leads to a contradiction (undamped oscillations 
are regarded as unstable). One cannot infer, however, 
that unstable oscillations are ruled out when g is less 
than the critical value at which static equilibrium be- 
comes unstable, because dynamic forces are not con 
sidered in the static investigation. Hence the critical 
value of g, above which the static equilibrium becomes 
unstable, is an upper bound to the critical dynamic pres- 
sure for panel flutter in the dynamic sense. Although 
the static considerations throw no light on the dynamics 
of panel flutter, this upper bound of critical dynamic 
pressure turns out to be alarmingly low under certain 
circumstances, and hence the result is significant. 

The same problem has been investigated by Isaacs,' 
Hayes,” and Miles.* 
change in length of the panel due to the axial thrust, 


However, they neglected the 


i.e., the strain in the mid-surface of the plate caused by 
the membrane stress is neglected. In this paper this 
change in length is not neglected. It will be shown that 
the membrane strain has far reaching effects on the 
stability of the plate, and entirely different conclusions 
are reached with regard to thermal buckling. Accord 
ing to the earlier analyses, a finite critical dynamic pres- 
sure for static stability exists regardless of the initial 
amplitude of the buckle. When the membrane strain 
is considered, however, this critical dynamic pressure 
will depend on the initial amplitude of the buckle. In 
fact, it tends to zero as the initial amplitude of the 
buckle approaches zero. 

The dynamics of panel flutter in the linear case of con 
stant membrane stress has been treated by Miles® 
(using approximate aerodynamic coefficients), and 
Shen‘ (using exact aerodynamic coefficients derived 
from the linearized aerodynamic equations). A relation 
of Shen’s solution with the rigid-support case of the 
present paper will be discussed later. 


(2) GENERAL EQUATIONS 


Consider a strip EF (see Fig. 1) of a two-dimensional 
panel as specified at the beginning of the Introduction. 
Let one edge of the plate be hinged, that is, it is free to 
rotate but is fixed in position, while the other edge is 
attached to a spring, with a spring constant a. When 
the spring-supported end is displaced by a distance A, 
the compressive force per unit length in the plate will be 


IT = II, + ad (1) 


where //, is the initial value of //. In the initial con- 
dition, H7) acts in the plate, the mid-surface of which 
can be represented by 


x 


ye 2 


m= 1 


. mr : 
a, sin (2) 


L 
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Under a lateral loading, the mid-surface becomes 


we 


RJ = 


- ici max (: 
m= 1 L 
These expressions satisfy the boundary conditions of 
hinged edges. 

When one side of the plate is exposed to a supersonic 
flow, the aerodynamic pressure p, acting on the plate is 


by = (2/V M — 1) q (dy/dx) (4) 


where g = (1/2)pl” is the dynamic pressure of the flow, 
U is the free-stream velocity, and / is the Mach num- 
ber. To this may be added a distributed load per unit 


area of the mid-surface, p: 


- -. MBX . 
bp=p Dd sin (5) 
m= 0 L 
The equations of equilibrium take a simple form if 
we assume 
lwl<L, Ilyi<«L, lanl <1, 
| dy . " 
lbn| <1, -} <1 (6) 
dx 


and that the thickness of the plate is much smaller than 
the radius of curvature of the plate. Then the well- 
known large deflection equation of von Karman, as 
generalized to initially curved plates, can be specialized 
into the following: 

d*yo 


d*y 
= = Hy - 
ax? dx? 


= d‘(y — Yo) 
D Eo apn 
dx* 


=—(p+pD,) (7 


where D = Et*®/12(1 — yw’) is the flexural rigidity of the 
plate, ¢ is the plate thickness, / is the Young’s modulus 
and p is the Poisson’s ratio. The difference in the mem- 
brane stress // from the initial value [/) is caused by the 
shortening of the centerline of the plate and the com- 
pression of the spring. Hence, assuming no compres- 
sive force acting in the direction of the edge AB or CD 


Et f | (2)" ()"] 
_ dx — 
oF, 0 dx dx 
Et 
S) 
L A ( 


we have 


H = Hy+ 


Eliminating A between Eqs. (1) and (8), and substitut- 
ing (2) and (3) for yo and y, we obtain 


roe) 


H = He +6 = Ym? (an? — bm?) ~— (9) 
m= il 
where 
= ——_— 10 
i + ae (10) 


Substituting Eqs. (2), (3), (4), (5), and (9) into Eqs. 
(7), an equation of equilibrium expressed in the form of 
Fourier series is obtained. On equating the coefficients 
of the corresponding terms in the right- and left-hand 
sides of this equation we obtain a set of simultaneous 
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equations 


‘Et . 
B vd >» n*(b,,2 — an”) | mb, + 
n 1 


tL 
[= nt Hoe) oe ee eee 
L* L 
2q , 2mn{l (—1)"*" Jd, 
Vi — 1. 4 (m + n)(m — n) 
eer 1 es a) 


where >>’ indicates that the term m = n must be ex- 


cluded. 
Introducing the dimensionless parameters 
An B. = V12(1 —- 
An = V12A1 — pw?) , 5 ¥ 
-* 2 b aL 
2t 
2L* q 
S = AyL?/x*D, QO = 
wiV M*?-—1D 
R V 12 eed l 
= *? ( — 2 ( 2 
12(] Ee) 2n'Dt 


Eqs. (11) become: 


B, (s - 
n=1 


n*B,? — B =. n°>r,,7 +m? — s) = 
mn = 


. ' R Q , 
\,(m? — S)- — kn - S&S D> 
m m? na 


2mn{1 — (—1)"7" JB, 


(m + n)(m — n) 


(m = 


The parameter X,, is proportional to the ratio of the ampli- 
tude of the mth harmonic of the initial arch to the plate 
thickness, B,, is proportional to the ratio of the amplitude 
of the mth harmonic of the deflected surface to the plate 
thickness, S is the ratio of the initial axial compression to 
the Euler buckling load of the plate, R is a dimensionless 
quantity specifying the additional loading, and Q is a di- 
menstonless quantity proportional to the dynamic pressure 
of the flow. Our problem is to determine the relation- 
ship between B,, and Q when R = 0. 

In the earlier works of Isaacs' and Hayes’ quoted 
above, the following differential equation is studied: 


d*y 2 dy 


dx? ' VM?—-1 


d‘y 
D— 


dx dx 


This equation can be derived from Eq. (7) by setting 
Hy = yo = 0. For a given plate, a pair of eigenvalues 
H and gq can be determined for which a nontrivial solu- 
tion exists. Without asking how J// is related to the 
amplitude of the deflection surface, the totality of such 
pairs of eigenvalues defines all the possible equilibrium 
configurations of a plate that is flat in the unstressed 
state. On setting y in the form e” and forming the 
characteristic equation according to the hinged edge 
conditions, Hayes gives a sketch of the eigenvalues simi- 
lar to that shown in Fig. 2. It is seen that the eigen- 
values form a series of closed loops. 
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Fic. 2. The first two loops of eigenvalues, relating the change 


of axial thrust to the dynamic pressure of flow (8 = 1, S = 1) 


However, if the supports were rigid, // varies with y 
and the more complicated situation must be clarified 
using Eqs. (13). 

Before discussing specific results, we shall give an 
alternative derivation of Eqs. (13) on the basis of energy 
consideration. Regard b,, as generalized coordinates, 
and let strain energy in a plate strip of unit width due to 
the axial compression // be U,, that due to the spring at 
the support be U,, and that due to bending by U,,. 


Then 


U, = H*L/2Et, U, = H?/2a (15) 


To calculate the strain energy of bending, it is important 
to note that yo, designated as the initial deflection, is not 
an unstressed state: it is the deflection under axial 
compression //). The unstressed plate will have a mid- 
surface y, given by 

D(yo” — yu") = —HAoyo 


Thus 


L 
U, = (1/2) [ D(y"” — yu")? dx = 
0 
L 
(1/2) f D(y” — yo" — (Ho/D)y)? dx (16) 
0 


Substituting Eqs. (2), (3), (9), and (10) into Eqs. (15) 
and (16), we obtain the total strain energy 


U, . U, + l s = 

Lf rkt < : : ae 
+ Ho 8 m? (dm? — bn2)e + 

ie} 2Et ' r 4 », 

AE rm HH? — , L 

2 V2, On — an) Tt pe Bang t 

| Pie m?a?L 

2 > an(bm — Am) - t (17) 
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The equation of equilibrium is then given by the La- 


grange’s equations 


[(o(U,, + U, + U,)]/0b,, = generalized force (18) 
which can be identified with Eqs. (13). 
(3) CasE 1. Fiat PLATE—INITIAL THRusT LEss 


THAN EULER BUCKLING LOAD (RIGID-SUPPORT, 8 = 1) 


If the compressive stress induced by the thermal ex- 
pansion is less than the fundamental Euler buckling 
load of the plate, the plate will remain flat at all values 
of dynamic pressure. In other words, the undeflected 
plate is the only statically stable configuration. 

To prove this statement, multiply Eqs. (13), with R 
= Oand # = 1, by m*B,, and sum all equations over m: 


Lim°B..?( Ln°B — Din*h,? + m* — s) = 


mr »,B»m(m*? — S) (19) 


The present case is characterized by the conditions 


a= 3 and Aj = 0, imme £2 Z....) 
which render Eq. (19) into 
(2mB,.2)* Lm 2(m? — S)B,”2 = 0 (21) 


Since every quantity on the left-hand side is positive, 
it is clear that the only possible solution is 


B, =0 


for all m (22) 
There exists no nontrivial solution, no matter how 
large Q is. 


(4) Case 2. FLAT PLATE—INITIAL THRUST EQUAL TO 
THE EULER LOAD (RIGID SUPPORT, 8 = 1) 


If the initial thrust is equal to the Euler load* the 
plate buckles into a sinusoidal arch, the amplitude of 
which depends on the amount of thermal expansion. 
To fix the idea let us assume A; > 0 so that the plate 


buckles into the airstream. In this case, 


S = l. Ay > 0, Ao = As =...=0 (23) 


When Q = 0, the only solution of Eq. (13), with R = 0, 
is actually the initial condition: 


B, = &, B, = O form # 1 (24) 


For increasing values of Q (necessarily positive if it is 
to have a physical significance) B, decreases and B,,’s 
become finite. From Eq. (19) one can deduce at once 
a necessary condition 

> m?B,.? < di? (25) 


m 


* Note that the initial thrust can never exceed the fundamental 
Euler load if the plate is to remain stable. 
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The interpretation of Eq. (25) is simple, for according to Hence ] 
Eq. (9) one obtains ati : 7 yas 
x = 3/2 + (3/2)V 1 — [(8/9)QO}? (31 sa 
S{(H/Hy)) — 1] = 2 — YS mB, (26) vari 
m x becomes imaginary when Q > 9/8, which is a critical J of (: 
Thus Eq. (25) implies that the aerodynamic load can Value for Q, the meaning of which will be discussed later. J is s¢ 
only cause an increase in the axial thrust in the plate. However, for the solution to have any physical mean It 
It is interesting to note that although the total lift ing, B, must be real. From Eqs. (28) and (30) one ob. J four 
acting on the plate is zero, for according to Eq. (4) the tains ishe 
total lift is proportional to y(L) — y(0), which vanishes ‘ a plet 
in our problem; the axial thrust does change because «ties 27x Q* I 
the elastic curve shortens under the asymmetric aero- mor 
dynamic pressure. Hence x must not exceed \,*._ This limitation on x has ior 
To solve the simultaneous Eqs. (13), with R = 0,and a very simple meaning. For when B, vanishes, all B, | "8! 
8 = 1, observe first that | B,, | decreases very rapidly as vanish and the plate becomes flat. Obviously no | @ is 
m increases: for according to Eq. (13) further increase in axial thrust is possible. The neces. | ™™ 
~y ae iti qua 
B,.| wr tihia® (27) sary condition ae 
. pee . , : x <A (32 : 
when m is large. The method of successive approxima- mate 
tion can therefore be applied. Let is, of course, also implied by Eqs. (26) and (28) and is | ™ a 
. ‘ as oe not due to the two-term approximation Eq. (30). eames 
aye ET ee Lm'B, ; -” For each Q the two fade a lead to two vt of B, tion 
(The sign of B, is determined by that of \; according to mad 
Then the equations to be solved in the present case are: Eqs. (24).) Having obtained B,, B, can be computed 
8 1G. 24 from Eqs. (30). The deflection of the arch as Q gradu- 
xB, + 3 QB: + 15 QBs + 35 QB +...=0 ally increases is thus determined. _ 
» 6 10 The variation of B, with Q is shown in Fig. 3 for the P . 
3 QB, + (x — 3)B, + 5 QB; + 7] OB; + case \; = 5. When Q = 0, B; = \; = 5 is indicated by | ™ ‘ 
= the point P,. As Q gradually increases, B, decreases 
; , (29) along the Curve P,P. At the point Po, Q reaches the 
= bad OB, + (x — 8)B; + bes OB, + critical value 9/8. For further increase of Q continuous 
15 21 change in the static equilibrium configuration is no = 
h.- OR +... = 0 longer possible. At P,: the two roots of x given by Eq. ee 
KF 7 ee ; 1.2 | 7 
‘or the existence of a nontrivial solution the determi- | 3 Pe | 
nant of the coefficients of B,, must vanish. This determi- i J 7) i 
nant A(Q, x) isa function of Q and x from which x can be LO a a \ ‘a i 
solved for a given value of Q. The order of the deter- I / \ / \ 
minant being infinity, there exists an infinite number of + / s ‘B 
roots x for each value of Q. Inasmuch as A(Q, x) speci- 0.8 18s } wn ! = eee ; ' 
fies the relation between the axial thrust and the dy- i! | / | : | 
namic pressure for static equilibrium of a panel whose r ; | an \ Q 
unstressed state is flat, it is equivalent to the equation Q 06 = ae / = ( 
determining the eigenvalues //, g of Eq. (14). The con- fi Moly | 
nection between // and y given by Eq. (8) determines | i / | | 
those eigensolutions that are physically realizable. A 0.4 ae Bo 4 4 ( 
further condition of stability then singles out the / | ! \| 
stable configurations. Pg | 7 A=5 | 
The characteristics of the solution are exhibited by | | a oe - ( 
the following approximate calculation. Neglecting all r | 
B,,’s for m > 3, Eqs. (29) become y Pe | P, 
| | 
xB, = —5QB, (x —3)Br= 2B, (30) . lh Uc nh hh 
3 3 Bn 
: Fic. 
The determinant of the coefficients is a Abii tance Gah te Spee sana ee Oe 
x? — 3x + (16/9)Q? = 0 a dan tan aa) Ne 
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(21) coalesce, hence upon reduction of Q the arch may 
pass on to the other branch P,0. The corresponding 
variation in B, is alsoshown. Using the third equation 
of (29), B; can be estimated and is plotted in Fig. 3. It 
isseen that B; is small in comparison with B, and By. 

In Fig. 4 is shown the variation of B,; with Q for 
four different values of \;. In the case A; = 1, B, van- 
jshes when Q = 1.08, at which the plate becomes com- 
pletely flat. 

The above calculation can be refined by retaining 
more equations in forming the characteristic equation 
When four equations are retained by 
., a fourth order equation for x and 


for x and Q. 
neglecting B;, Be, . . 
( is obtained. The roots of this equation are deter- 
mined numerically and plotted in Fig. 2. (If more 
equations were considered, further larger loops would be 
added farther to the right.) 
Note that the maximum value of Q 


The first loop is approxi- 
mated by Eq. (32). 
in the first loop is approximately 9/8. Since this peak 
value of Q is of some importance, a further approxima- 
tion by taking six equations (neglecting B;, Bs, . . .) is 
made, which gives 


Oer = 1.172 (33) 
when A; > 1.255. If Ay < 1.255, the critical value of Q 
at which the plate becomes flat depends on the value of 


\; and 
(34) 


Ov, — 0 when \; ~ 0 


(5) CONDITIONS OF INSTABILITY 


The question of stability may be approached from the 


point of view of potential energy. If the potential 
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of the deflection curve with the dynamic pressure of the flow, 
lor several values of initial amplitude (A,). Dotted lines indi- 
cate unstable configurations. (Four term approximation, 
B=i1,S = 1). 





PANEL 


IN SUPERSONIC FLOW 561 
energy of an equilibrium configuration is a relative 
minimum with respect to neighboring configurations. 
the equilibrium is stable; otherwise it is unstable. 

We must point out that the theorem of minimum po- 
tential energy as presented in most elasticity books does 
not apply to our problem, even though our plate re- 
mains elastic. Indeed, if the truth of the theorem is 
not qualified by auxiliary conditions, it would have im- 
plied that all equilibrium states are stable, to which an 
obvious contradiction can be given by an example; 
namely, a straight column subject to an axial thrust 
greater than the Euler buckling load. For such a col- 
umn a slight external disturbance will cause buckling to 
take place. The undeflected configuration, though in 
equilibrium, is unstable. 

The proof of the minimum potential energy theorem, 
as presented in reference 5, assumes that the elastic 
strain can be computed from the displacements by the 
), where 1, Uo, U3 


linear relation e;; = (1/2)(u;, ; + uj. 9 
are the components of displacement with respect to a 
fixed Cartesian coordinates system x, Xs, x3. This is 
not the case in the derivation of von Karman’s equation 
from which Eq. (7) is obtained. For example, the 


strain component e,; in the plate is 


Ou 4 1 (S*) 1 (Se) 
eé,- = om» 
7 ox * 2 lox 2\or/ 7 


where y and y, are the vertical deflections of the mid- 


aXX(y — y,) 
Ox? 


surface of the plate in the final and unstressed states 
respectively, z is the distance from the mid-surface, and 
u is the displacement of a point on the mid-surface in 
the x-direction. The nonlinear strain-displacement re- 
lationship affects the variational process in such a way 
that the theorem of virtual work still holds, while the 
potential energy may not necessarily be a minimum. 
The general statement that the minimum potential 
energy theorem is “‘capable of general application in 
that it is true whatever the law connecting load and 
deformation may be,’’ as can be found in some publica- 
tions, is quite erroneous. 

A simple criterion of stability is obtained by consider- 
ing the force that is required to balance a system that is 
disturbed from an equilibrium position. If this force 
tends to push the structure back to the original position, 
the equilibrium is stable. Otherwise it is unstable. Al- 
ternatively one may compute the work done by the ex- 
ternal force that is required to balance the disturbed 


structure when the structure is disturbed from the 
equilibrium position. If this work is positive, the 
equilibrium is stable. Otherwise it is unstable. Note 


that if the work referred to above is positive for all 
possible disturbances which satisfy the boundary con- 
ditions where displacements are specified, the potential 
energy of the system is a minimum. 

Applying this concept to our plate, let B,, Q, and x 
describe an equilibrium configuration. Let B,, be given 
an arbitrary small disturbance 6B,,. Then Eqs. (13) 
give the external force (described by 6Rk,,) that is re- 
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required to balance the structure in the disturbed 


position : 


(Bn + 6Bm) X 
E > n(B, + 6B,)? — B 
1 


n 


>. n*r,7 + m*? — s| _ 
i = 1 


) 
An(m? — S) + . 


m* 
2mn[{1 — (—1)" * "](B, + 6B,) oR 
(m + n)(m — n) m* 


n l 


’ 


(35) 


Subtracting from this the original equation in which 
6B, = 6R = O, and neglecting (68,,)* in comparison 


with 6B,,, we obtain 


6R 
ad kn = 26B, >, n°B,6B, + 


m* Big 
éB.. 18 = n°B,” — B >» nv," m— S|+ 
n 1 n= 1 
y. 2. See lt = (= {)" * * 
e : 5B,, (36) 
m*,=1 (m+n)(m — n) 


The external pressure distribution is 


2r'Dt marx 


p=7 = 6bR > k, sin 
V121—y2)L' 4 cn) ; 
The work done by the external pressure through the 
disturbance, on a strip of unit width, is 


1 L 
~ f p(x)by(x) dx = 
< Jo 


w*Di? 
12(1 — p*)L? 


W = 


6R > knbBm (37) 


m= 1 


From Eq. (36), we obtain 


W = {wDe/[12(1 — u2)L*)} 12 m(B,,)? X 
[3 hs n*B,? — B u n*r,2 + m?>—-— S+ 23m°B,,? | + 
= p ig 4Bm°n°B,,B OB 5B »\ (38) 


m n 


This is a quadratic form in 6B,,, (m = ..) and can 


| ie ae 
be written as 


{rtDt2/(1211 — w2)L*]} YS an,bB, 6B, (39) 


mn = 1 


W = 


The coefficient of Q vanishes in Eq. (38), a fact that can 
be explained by Eq. (4): a symmetric deflection (with 
respect to the center of the panel) of the plate induces 
an antisymmetric aerodynamic pressure distribution, 
and vice versa, so that the integral over the span of the 
product of pressure times deflection vanishes. 

The condition that the plate be stable is that the 
quadratic form IV be positive definite. A set of neces- 
sary and sufficient conditions for this is 


ay, > 0, lau ay| > 0, Qy; Ay Ay3| > 0, Rae 
ay Ax Az, Az2 = a3} 
(131 39 (133 | (40) 
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where 

a m? X 

(s >> n°B,? —B > nr,2 + m2 — S + 28m°B,?) 
n= 1 n l / 


AGnn = 48m°n7B., Bn, (m # n) (4] 


Apply these results to the case studied in the last section 
in which 8 = 1, S = 1, A; + 0, and A, = Ofor m #0 
we obtain the conditions for the plate to be stable: 


— x + 2B,? > 0, 


|—y + 2B,2 16B,B, +0 
16B, By 4(—x + 3 + SB,’) 
|_y + 2B, 16B,B, 36B,Bs $0 
1GB,B,  4(-x + 3+ S8B.2) 144B,B; 
36B,B; —-144BB; (—-y +8 +4 
18B;2) 
(42 


etc. 
If any of these inequalities are violated, the plate be 
comes unstable. 

It is now clear that the stability of the plate con- 
sidered in the last section depends on the initial ampli- 
tude \;, because for each pair of eigenvalues x, Q, the 
ratios of B,, Bo, ...can be solved from Eqs. (29) so that 
a unique deflection mode is obtained, but the amplitude 
of the deflection mode must be solved from Eq. (28 
which depends on )). 

When the stability criteria, Eq. (42), are applied to the 


equilibrium configurations represented in Fig. 4, we 


find that instability sets in long before the value Q,, = 





TABLE 1 


Eigenvalue, Mode Shape and Stability for the First Loop 












































- | | 
; oe}; x | r, i | B, | 33 | % 
}0.3 | 0. 05461 | 1.0 0.964 | -0.065 i 0.001 -0.001 Stable 
2.5 2.466 | -0.167 | 0.003 | -0.003 | Stable 
5.0 4.949 | -0.335 | 0.007 | -0.006 | Stable 
10.0 9.906 | -0.671 | 0.013 | -0.013 | Unstable 
— +} ____ + 
lo. 35 | 0.07479 | 1.0 0.950 | -0.076 | 0.002 | -0.002 | Stable 
| | 2.5 2.454 | -0.195 | 0.004 | -0.004 | Stable 
; 5.0 4.930 | -0.392 0.009 -0.008 | Unstable 
L 10.0 9.872 -0.785 0.018 -0.016 | Unstable 
| T a _ 
0.5 | 0.1565 | 1.0 0.894 | -0.104 | 0.003 | -0.002 | Stable 
2.5 2.404 | -0.809 | 0.009 | -0.005 | Stable 
} 5 0 4.854 | -0.565 0.019 -0.011 Unstable 
| |} 10.0 S738 | «3.333 0.037 -0.023 Unstable 
}+——_+ . j+—— + $$ $$ 
}0.6 | 0.2304 | 1.0 0.843 | -0.121 | 0.005 | -0.002 | Stable 
| Z.5 2.359 -0. 337 0.014 -0.006 | Unstable 
5.0 | 4.784 | -0.684 | 0.027 | -0.013 | Unstable 
| 10.0 | 9.602 | -1.372 | 0.055 | -0.025 | Unstable 
—s i | 
0.9 | 0.5841 1.0 | 0.581 | -0.140 | 0.004 | -9.003 | Unstable 
| 2.5 | 2.143 | -0.518 | 0.015 | -0.01 Unstable 
| | se | 4.448 | -1.074 | 0.031 | -0.020 | Unstable 
10.0 8.975 | -2.168 | 0.063 | -0.041 | Unstable 
ss | | 
}1.0 | 0.7748 | 1.0 0.411 | -0.118 | 0.009 | -0.002 | Unstable 
2s 2.024 | -0. 583 | 0.042 | -0.010 | Unstable 
5.0 4.258 | -1.226 | 0.088 | -0.023 | Unstable 
| 10.0 | 8.618 | ~2. 482 | 0.179 | 0.046, | Unstable 
r + 4 
}1l.1 | 1.0545 1.0 | imag | | | 
| 2.5 | 1.850 | -0.660 | 0.054 | -0.012 | Unstable 
| 5.0 3.972 | -1.142 0.117 | -0.025 | Unstable 
10.0 8.073 | -2.882 | 0.237 | -0.050 | Unstable 
90.0 | 73-04 -26.07 | 2.147 | -0.456 | Unstable 
= 4 
fi. [ 2-10 | 1.0 imag | | | 
2.5 1.161 | -0.828 | 0.081 | -0.009 | Unstable 
| | oe 2. 108 -1.944 | 0.190 | 2023 | Unstable 
9 .641 | -4.020 3 - ble 
| | 0-293 | -9.047 | Unsta 
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Eigenvalue and Mode Shape for the Secor 





TABLE 2 


id Loop 


Value of B) for A, = 


(all unstable) 





PANEL 








1 3 4 1 : tie ts 
2.5 5.0 10.0 

— + + + + + - — 
1.0 7.923 -2.608 11.6 0.71 -* 0.117 0.272 
14.95 -0.10 1. 50 -13.8 : 0.057 0. 167 
2.0 7. 746 30 3.12 0.39 - 0.42 0.97 
14.79 -0. 178) 1.43 -6.49 : 0.12 | 0.35 
— La Sa Ee ———+—— + + Ee 
4.0 7 692 -0. 59 | 1.14 -0. 30 - 1.08 2.43 
14.11 -0.28 | 1.20 -2.61 ° 0. 30 0. 84 
— $$ +—_——— — 
5.0 8.054 -0. 47 | 0.95 -0.34 - 1.20 2.78 
13. 51 -0. 30 1.08 -1.78 - 0.43 1.18 

———+ — +— ae | cae | ~ 4 4 — 
6.0 8. 84 -0. 38 0. 86 -0.43 ° 1.20 2.85 
12. 56 -0. 31 0.97 -1.18 ° 0.62 1.65 

_ T — T —i — + — + ———+>—-———_ 
6.5 9.62 -0.35 0.85 -0. 52 - - 1.12 2.71 

1.71 0. 32 0.91 -0.90 ° ° 0.78 2.01 | 


-— a4 


No root < 24 





* A formal solution gives an imaginary B) 





1.172 is reached. In Fig. 4 the stable configurations are 
drawn with solid lines, while unstable configurations 
are shown dotted. Furthermore, for these \, values all 
equilibrium configurations corresponding to the eigen- 
values on the second loop in Fig. 2 are unstable. Hence 
when Q exceeds the value given by the points joining 
the solid and dotted curves in Fig. 4, indicated hereby 
as Q,, instability sets in. 

Tables | and 2 show the eigenvalues Q and x; the 
equilibrium configurations specified by the Fourier co- 
eficients B,, B., B;, and By, (B,,’s for m > 4 are neg- 
lected); and the stability of the equilibrium configura- 
tions. The values of x corresponding to each Q in Table 
| are the smallest eigenvalues, except for Q = 1.1 where 
both the first (x = 1.0545) and the second (x = 2.10) 
eigenvalues are given. All equilibrium configurations 
corresponding to the second and higher eigenvalues of 
x (for each specified Q) are unstable. As to the equi- 
librium configurations, in general, the first harmonic 
B,) predominates in the first mode (the smallest eigen- 
value x), the second harmonic (B2) predominates in the 
second mode (the second eigenvalue x) and so on. 

The critical value Q.,, above which statically stable 
equilibrium does not exist, is shown in Fig. 5, as a func- 
tion of A). 


(6) INITIALLY CURVED PLATE 


In Sections 3 and 4 the thermal gradient in the plate 
thickness is neglected. With thermal gradient the 
plate will be bent even though the axial thrust is less 
than the Euler load. The exact initial curvature and 
axial thrust depend on the conditions of heat flow, the 
calculations of which will not be discussed here. It is 
assumed that the initial curvature and axial thrust are 
known. 

As an example consider the case 


B8=1 
(43) 


1, ¥ 0, Ww =A=...=0 and 
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for which the governing equations reduce to Eqs. (29) 


if S = 1. If S < 1, the variation of B,, with Q can be 
traced according to Eqs. (13). Again as an approxima- 


tion let us neglect all B,, for m > 3, so that the funda 


mental equations become 
i1-S 
S)B. = 0 14 


(x — 1+ S)B, + (8/3)QB, = —r 
(2/3)0B, + (x — 44 


where x is defined by Eq. (28). Solving Eq. (44) for 

B,, B,; substituting into Eq. (28) and writing 
16 Q? is 

n = 2 | »~) 

94-S-, 

one obtains: 

Ay"(1 — SS)? 
X= A? = : (1 + : ) (46) 
(x — 1+ 5 — »)* i—-S— x 


For given values of \;, S, and x, Eq. (46) can be solved 
for », which in turn can be converted into Q? by Eq. 
(45). Thus the complete dynamic pressure vs. deflec- 
tion curve can be traced by gradually increasing x, 
finding the corresponding yn and Q, and finally B, and 
B,. Such a calculation for the case \; = 5 yields the re- 
sults shown in Fig. 6. As in the previous section B, 
decreases initially with increasing Q, and a maximum 
value of Q exists in each case. This maximum, denoted 
by Q.,, is plotted in Fig. 7 as a function of .S. 
that Q., increases almost linearly as the initial thrust de- 
0, i.e., without initial thrust, 


It is seen 


creases, and that at S = 
Q., is of the order 1.7. In Fig. 6 are also shown a point 
for \; = 10, S = 0; and a point at S = | given by the 
results of Section 4; thus indicating the effect of \,; on 
Q., and also the effect of the terms neglected in the cal- 
culation. 

When \, is small and S < 1, the plate tends to behave 
like an For example when S = 0 
and \, 


unbuckled flat plate. 
= |, the axial thrust in the plate amounts ex- 














ll 


The variation of Q., with \; (8 = 1, S 


Fic. 5 
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B, 3.0 4.0 


Variations of axial thrust and dynamic pressure with 
deflection, for an initially curved arch with \; = 5. 


Fic. 6 


actly equal to the Euler load when the plate is made flat. 
In this case there is no critical dynamic pressure: the 
plate becomes flatter as Q increases, and becomes a 
plane as Q becomes sufficiently large. 

The stability of the equilibrium configurations can be 
screened again by means of Eqs. (40) and (41). The re- 
sults are indicated also in Figs. 6 and 7, where unstable 
The 
actual critical dynamic pressure, given by Q,,, is seen to 
be much lower than Q,,.. Q., gives the dynamic pressure 
above which statically stable equilibrium does not exist. 


configurations are identified with dotted lines. 


(7) CONCLUDING REMARKS 


The behavior of the plate at very small \; is such that 
Q., > 0if S=1, while Q., > © if S< 1. This peculiar re- 
When 
dynamic considerations are made this singular behavior 
of the plate as S— | and \,; ~ 0 is removed. Since the 
effect of the rigidity of the supports tends to vanish as 
\i > 0, the analysis of the dynamic stability of the plate 
can be made disregarding the restraint imposed by the 


sult is derived from static considerations alone. 


fixed hinges. The analysis by Shen‘ is applicable to 


| of ff | 
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Fic. 7. Critical dynamic pressure as a function of initial axial 


thrust, for initially curved plates (6 = 1). 


plates with constant membrane stress; inasmuch as the 
effect of membrane stress consists of only a change in 
the natural vibration frequencies of the plate. In ap- 
plying Shen’s* analysis, let us notice that when S > 1, 
\, > 0, the natural vibration frequencies of the plate for 
infinitesimal amplitude of vibration are w, = 0 for the 
first mode and 


we = W12(x2/L2VEI/o (47 


for the second mode, where J and o are the moment of 
inertia and mass per unit length of the plate, respec- 
tively. An investigation of Shen’s results shows that 
under the stated circumstances the critical dynamic 
pressure of flutter tends to zero. Thus the statically 
stable configuration when ), is small and S is slightly 
less than | is dynamically unstable, just as when S 1s 
equal to 1. 

When Q > Q,,, static equilibrium being unstable, 
motion ensues. This motion cannot be a stable oscilla- 
tion (undamped oscillating being regarded as unstable). 
Hence Q,, is an upper bound to the critical dynamic 


pressure for panel flutter. This upper bound is shown 
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STABILITY OF CURVED 


tending to zero as d; tends to zero. The practical im- 
portance of this limiting case lies in its relation to ther 
mal buckling: in transient heating the stage A; = O is 
If the 
heating is slow the plate may stay long enough in this 
Hence one must con- 


always reached before ); passes to a finite value. 


danger period to cause flutter. 
clude that no thermal buckling can be permitted for any 
fat panel of a supersonic aircraft, unless there is reason 
to believe that the heating is such that the buckling 
amplitude jumps at once from zero to a sufficiently large 
magnitude. 

In the above only the “‘rigid-support’’ case is con- 
sidered. In practice both the “‘rigid-support’’ case and 
the ‘‘perfectly-flexible-support’’ case are important; 
the former being a limiting case of a panel attached to 
heavy frames; and the latter, a limiting case to dia- 
phragm type of supports. In the latter case the mem- 
brane stress (axial thrust) may be regarded as a con- 
stant independent of the amplitude of the deflection. 
The equilibrium configuration of such a bent, flexibly- 
supported plate can be easily computed. 
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(Continued from page 524) 
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Stability of Inelastic Systems 
(Continued from page 548) 


Pearson, C. E., Bifurcation Criterion and Plastic Buckling of 
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17, pp. 417-424, 1950. 
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The more general case of an elastically supported 
panel can be solved by the same method given above, 
by giving 8 its proper value. The effect of finite angle 
of attack can be investigated by incorporating proper 
Rk,, terms. 
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On the Excitation of Pure Natural Modes 


Ernst Schultze 
Mathematician, Flug- & Fahrzeugwerke A.G., Altenrhein, 
Switzerland 
March 8, 1954 
| Regeaes AND WRISLEY have new experimental 
method of exciting pure natural modes of complex struc 
This 
method is based on the obviously experimentally stated fact, 
that for the excitation of a pure natural mode, all the exciting 
forces must be in phase and all the force amplitudes must be pro 


presented a 


tures, such as airplanes, by a system of many shakers.! 


portional to the product of the amplitude of that natural mode 
and the mass distribution, in formula, F/(s-m), independent of 
individual shaker. Lewis and Wrisley attempt to explain this 
fact by aid of structural damping. 

Considering the great value of the method of Lewis and Wris 
ley, which permits a very quick adjustment of pure natural 
In the fol 


lowing, such a proof is given for a nearly plane wing on the basis 


modes, an exact proof of the above fact is desirable 


of the theory of integral equations with symmetrical kernel 


SYMBOLS 


F, Q, points on the wing 
2(P) complex amplitude of P 
m(Q) mass distribution 
dQ = dx dy wing surface element 
m(Q) dQ mass of element dQ 
6(P, Q) elastic influence function deformation at P when 

force at Q is unit) 
i time 
v circular frequency 
2 coefficient of structural damping 
0; point of attachment of shaker 7 (j 1.2, n) 
Fj force amplitude of shaker j 
Vi, v2, 3, natural frequencies of wing 
2i(P), 22(P), natural modes of wing 
ny parameter of integral equation 

/ 
¢(P) 2(P) V m(P) 
K(P, Q) 5(P, 0) V m(P)m(Q) 
k(Q;) Fj/V m(Q;) (some abbreviations 
n 
F(P)  — ie) D> KP, 0,0) 
j=l 


DISCUSSION 
The vertical motion of the point Q being described by 
2(Q)e™ 
the inertia loading at the point Q is given by 
v?m(Q)s(Q)e'™dO 
This loading gives at P the deformation 


dz(P)e’ = v5(P, Q)m(Q)2(Q) dQe™” 


or by integration over the wing surface 


oP) = v? f° &P, O)m(Q)x(Q) dO 

This integral equation holds for the free oscillations of th 
wing without damping. With structural damping, the influence 
function may simply be multiplied by the complex factor (1 - 


ig), thus giving the integral equation: 


o(P) = v2(1 — 7g) - 6(P, Q)m(Q)s(Q) dQ 9 
If the wing is excited by a certain number » of shakers in th 
, n), the additional terms 


Ry 


points Q; (j = 


n 


S* (P.O 
2 0; 


1 


- Fe 
give the motion of the point P due to the excitation of the shakers 
Adding these terms we obtain the following integral equation for 
the forced vibration: 


. 
v? ff &P, Q)m(Q)“Q) dQ 4+ 


n . 
> WP, Q,)F, 
j=1 

Transforming now to a symmetrical kernel K(?P, Q), we intro 


F For the free 


duce the function ¢(/) defined in the list of symbols. 
vibration without or with damping, there results a homogeneous 


integral equation: 
? id ? ‘ 
oP) —-r» Sf K(P, Q)¢(Q) dQ = 0 
and for the forced vibration an inhomogeneous equation 


—r f K(P, Q)¢(Q) dQ = F(P) 


» 
¢g(/ 


The kernel A(P, Q) is not only symmetrical but also positive 
definite. 

In relation to this sort of integral equation there exists 
elaborated theory,? which vields among others the following 
facts: 

(1) All characteristic values (supposed simple), Aj, Az, As, 
of the homogeneous Eq. (4) are real and positive. 

(2) All characteristic functions, ¢i(P), g(P), ¢3(P), , ol 
the homogeneous equation are real, and, if normed, the following 
relations are valid: 

Pa ¢i(P)edtP) dP = : <f + : 

(3) If X + AX, it is possible to express the solution of the 1 
homogeneous equation in terms of the solutions of the homo- 


geneous equation: 


; fi, 
oP) = F(P) +2 >» gi(P) y 
ane 
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where 
f. = f° F(Q)exQ) dQ 


In the case of the free oscillating wing, we conclude: 

From (1), that the wing can only perform undamped (by ab- 
sence of structural damping) or damped oscillations with the 
frequencies ve = V Aj, k= 1,2, 3, 


4 
cies 


(natural frequen- 


From (2), that the natural mode z;(/) which corresponds to 
the natural frequency v, is a real function of the point P, which 
means that all points of a wing oscillating in a natural mode have 
the same phase. 

In the case of the forced vibration, we conclude from (3) after 
transformation, that the complex amplitude of the vertical motion 


of the point P is given by 


5 ay/(vy? — v®? + igv,?) |x(P (8) 


where 
an = > 2 k(Q;)¢Q;) = 
] ' 


This formula shows that in general all the natural modes of the 
wing are excited, and every natural mode appears with another 
phase angle. In the special case of the exciting frequency coin- 
ciding with one of the natural frequencies, the phase shift between 
the exciting forces and this natural mode is 7/2 (‘‘in phase’’ con- 
dition), the other natural modes have, however, other phase 


ingles 


Coming back now to the excitation of only one natural mode 
mentioned in the introduction, all the a,;'s in Eq. (8) except one, 


saV dm Must be zero: 


nm 
a, = >, k(Q/)e(Q;))=90 ifk=m (9) 
j=l 


But, if we have many shakers, the sum may be thought of as an 
ipproximation to an integral in the case, where the exciting forces 
were continually distributed over the wing: 


. k(Q;)¢%(Q;) ~ Dy (dk/dQ) ¢i(Q) dQ 
=] 


j 


The condition (9) is therefore approximately equivalent to the 


ondition 

S (dk/dQ) ¢(Q) dQ = 0 ifk +m 
vhich states that dk/dQ must be orthogonal to all ¢(Q) except 
gn(Q), thus dk/dQ proportional to ¢gm(Q) or approximately, if 
there is a finite number of shakers instead of the continually 
distributed exciting forces: 
[em(Q = F;/{zm(Q;)m(Q;)] = independent of 7 
This is exactly the claim of Lewis and Wrisley for exciting a 
pure natural mode. But now we see that the presence of struc- 
tural damping is not necessary to come to this conclusion, be- 


iuse Eq. (9) is independent of structural damping 
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Effects of Afterbody Configuration on the 
Pressure Distribution over the Nose of 
Two-Dimensional Wedges 


R. W. Truitt and F. W. Martin 
Virginia Polytechnic Institute, Blacksburg, Va 
March 3, 1954 


HE THEORETICAL incompressible pressure distributions over 
the front portion of two-dimensional wedges, with different 


afterbody configuration, have been calculated by the Schwarz 


Christoffel transformation. ! The cases considered here, at 
zero angle of attack, are: (1) the symmetric diamond wedge 
and (2) the single wedge with straight afterbody. See Fig. | 


The theoretical pressure distributions for these two cases are 
shown in Fig. 2 for a wedge semiangle, 6, of 5 It is interesting 
to note that there is a considerable difference in the pressure dis 
tributions for the two cases shown. In particular it is found 
that, for case (1), the Cp, = 0 point is at x/c = 0.32, whereas for 
case (2), the Cp, = 0 point occurs at x/c = 0.52 

In order to compare these theoretical results with experiment 
low-speed (Wo 0.1) pressure distributions were obtained for 


the two afterbody configurations. The experimental results, 


shown in Fig. 2, are seen to be in agreement with the theoretical 
results. The theoretical and experimental evidence shows the 
effect of these changes in afterbody configuration, namely, in 
creasing the change of flow direction at the shoulder causes the 
front-wedge pressure distribution to become increasingly more 


negative 
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For given free-stream boundary conditions, then, the incom- 
pressible potential flow over the front wedge is different for the 
two afterbody configurations. It follows, therefore, that the 
application of compressibility corrections to cases (1) and (2), for 
a given subcritical free-stream Mach Number, should yield en- 
tirely different pressures for the two cases. That is, in view of 
the results shown in Fig. 2, one should not expect the same sub- 
critical pressure distribution for the two cases even though the 
wedge angle is the same. 

Though no attempt is made, in the present note, to generalize 
these results, it seems evident that any change in afterbody con- 
figuration should result in changes in the incompressible front- 
wedge pressure distribution. Therefore, indiscriminate com- 
parison of front-wedge pressure distributions, based only on wedge 
opening angle and ignoring afterbody configuration, should be 
avoided, at least in the subcritical range. 

It seems reasonable to anticipate that the evolution of the 
pressure distribution from low speeds to critical and supercritical 
speeds, i.e., in the transonic range, should reflect the effect of the 


different afterbody configuration. 
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Comments on ‘‘The Effect of Boundary Layer 
on Sonic Flow Through an Abrupt Cross- 
Sectional Area Change’’! 


Helmut H. Korst 

Professor of Mechanical Engineering, University of Illinois, 
Urbana, Ill. 

March 12, 1954 


AS A RESULT OF recent work on base’ pressures in transonic and 
supersonic flow it is now possible to evaluate Wick’s experi- 
mental results on the basis of a purely theoretical analysis.? 
This analysis utilizes a two-dimensional flow model which con- 
siders the interaction between dissipative flow regions and the 
adjacent free stream in greater detail than discussed by Crocco 
and Lees.* 

Specifically, the physical principle of conservation of mass in 
the dead-water region has to be satisfied by formulating an ‘‘es- 
cape concept” applied to a particular streamline in the jet bound 
ary. If no mass is added to the wake externally, the 
concept” requires that the mechanical energy level of the stream 


““ 


escape 


line separating the approaching flow from the dead water is such 
that only the fluid of the approaching flow will penetrate into the 
high-pressure region created by the trailing shock. Under these 
conditions a unique and stable solution is obtained for the base 
pressure problem. A detailed analysis of the flow configuration 
within compressible jet boundaries, including information on the 
mechanical energy distribution, is therefore necessary. 

Such a study has been undertaken with special consideration 
being given to the influence of various velocity distributions at 
the beginning of the mixing region.4~® The analysis, which will 
hold in its present form as long as the jet boundary thickness 
remains small as compared to a characteristic dimension of the 
free jet, leads to an asymptotic velocity distribution which is 
gradually approached with increasing length of the mixing re- 
gion regardless of the initial distribution. This asymptotic solu- 
tion of the mixing problem representing the fully developed jet 
boundary corresponds to a lowest limiting base pressure for a given 
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Fic. 1. Comparison of experimental results with theory for 


external flow. 


Mach Number of the approaching free stream in any given con- 
figuration of the two-dimensional flow boundaries. Calculation 
of this limiting base pressure is possible without any empirical 
information. 

For external flow past sudden expansions, such a limiting be 
havior is reported by Eggink’ and is also indicated by an anal 
ysis of data compiled by D. R. Chapman, et al.8 Chapman's 
experimental results, obtained for thin approaching boundary 
layers, and Eggink’s ‘‘asymptotic values’’ are compared with the 
purely analytical solution for a fully developed jet boundary in 
Fig. 1. 

The theoretical solution for a fully developed jet boundary, 
as applied to the internal flow problem studied by Wick, is shown 
to be in agreement with his experimental data for shortest ap- 
proach length. See Fig. 2 

According to the two-dimensional jet mixing theory ther 
should be an increase in base pressure as the flow profiles in the 
mixing region are increasingly affected by the initial boundary 
layer. This result is supported by experimental investigations 
with external flow under more nearly two-dimensional conditions 
than those prevailing in Wick’s experiments and it must be con- 
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theoretical solution. 





clud 
repr’ 
the 

Wicl 
lime 
the 1 
wall 
theo 
pres 
rang 
initic 


three 


W 
pp. 67 
2Ci 
major 
iran 
Sectio 
Divisi 
Cr 
Dissif 
Scienc 
‘Cl 
derai 
of Illi 
5 Ch 
fInit 
Applic 
6 Ke 
Vol. 7 
Eg 
unne 


§ Ch 


The f 
1) in 


tions, 


al | 
0 


ory for 


Cll con- 
ulation 
1pirical 


ing be- 
1 anal 
pman’s 
indary 
ith the 
lary in 


ndary, 
shown 


‘st ap- 


there 
in the 
indary 
ations 
litions 


e con- 


4.0 


flow 
with 





READERS’ 


cluded that the adverse trend exhibited in Wick’s paper is not 
representative of the influence of the initial boundary layer on 
the two-dimensional base pressure problem. The decrease of 
Wick’s base pressure values below the theoretical limit for two- 
jimensional flow apparently is due to a decrease in strength of 
the trailing shock, caused mainly by the thickening of the side 
wall boundary layer. This latter effect is well known and a 
theoretical analysis shows that it will lead to a decrease in base 
pressure. . Since all of Wick’s experiments were most likely in the 
range of a fully developed, thin jet boundary, the influence of the 
initial boundary layer reported by him seems to be essentially a 
three-dimensional effect. 

It has probably come to the author's attention that Eqs. (2 


and (3) are erroneously printed. They should read: 


dR/R = (1/2)[(k + 1)/(k — 1)] tan ede —(1/2) ctnede (2 


R,/Re = (sin e/sin €,) /? (cos &/cos e )* (5 
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On Use of the Reverse Flow Theorem in the 
Problem of Aeroelastic Deflections 


H. M. Voss 
nstructor, Department of Aeronautical Engineering, Massachusetts 
Institute of Technology, Cambridge, Mass 


March 14, 1954 


HE INTEGRAL EQUATION for the deflection of any point on an 


elastic wing may be written 


Ax, vy, t) = JS Pi C(x, y; & n)gl(é, n, t) dE dn + 
S 


ie x C(x, y; & n) plé,n, t) dE dn (1 
s 
where 
(x, y, t) = deflection at any point 
C(x, y; & 7) = flexibility influence coefficient 
S = planform area 
g(é, n, t = inertia loading plus any impressed loading 
p(t, n, t = pressure loading due to motion, i.e., 
4 Lo Oz : 
w(t, n, t) = (&, 9, t) + (&, 7, t) (2 
, ot dé 


The fact that p is itself a complicated function of zs leaves Eq 
1) in a form unsuited to solution without simplifying assump- 


tions, e.g., strip theory aerodynamics or assumed modes 
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However, use of the reverse flow theorem! ? permits the re 


writing of Eq. (1) for steady or simple harmonic motion, in the 


form, 


(x, ¥, t) = / / C(x, y; & nelé, 
S 
. * (a , 
(&, 9, #) + (f, 7, t)} x 
ol of 


V(x, 9; &, nidE dn (3) 


n, t) dé dn + 


where V(x, y; & 7; t) = V(x, y; &, n)e'’ is the pressure distribu 
tion on the planform in the reversed flow due to a non-dimensional 


downwash distribution 
w(t, 9, t) = — é (4 
with w the circular frequency. It is to be noted that V is inde 
pendent of the deflection, z. 
From Eq. (3) V may be interpreted as a downwash-equivalent 
load function, giving the equivalent load at x, y per unit down- 
(3), is 


wash at £, 7». The reverse flow theorem, and hence Eq 


valid for any Mach Number, within the limits of the usual linear- 


ized theory 
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Slender Bodies of Minimum Wave Drag* 


W. T. Lord and E. Eminton 
Royal Aircraft Establishment, Farnborough, Hants, England 


March 10, 1954 


B be, WAVE DRAG OF A BODY slender in the sense that the slopes 
of its surface are everywhere small compared with unity may 
be calculated by linearized theory for a wide range of Mach 
Number. If / is the length of the body and S(x) is the distri 
bution of cross-sectional area for 0 < x < /, then for a body such 
that S’(x) is continuous for 0 < x < /, S’(O) = Oand S‘(/) = 0, 
the zero-lift wave drag D is given by 


D1 i hoa ! 
— S"(x) S”’(E) log dx dé, 
q 2r Jo x—¢ 
where g is the kinetic pressure. If the areas at the nose x = O and 


the base x = / are not both zero, and all the air in the stream tub 
defined by the nose area is assumed to pass through the body, the 
value of D given by this expression must be interpreted as the 
wave drag of the external surface between nose and base; the 
cylindrical for —© «< 


body may be considered to be 
x<Oandi <x 


on the cross-sectional area distribution, 


The wave drag coefficient depends only 
ind is independent of 
both the cross-sectional shape and the Mach Number 

Often only a few parameters defining the shape of a body: need 
to be specified, without the detailed shape being pre-determined, 
and it is possible to find the area distribution which minimizes 
the wave drag under these conditions. The authors listed in the 
references have derived the optimum area distributions for a 
variety of conditions. Now, paradoxically, the more restrictive 
the conditions imposed, the more general the resulting optimum 


* Acknowledgment is made to the Chief Scientist, British Ministry of 


Supply, for permission to publish this note 
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distribution, since the optimum distributions with fewer condi 
tions may be deduced by minimizing the wave drag with respect 
to those parameters which are no longer fixed. The optimum dis- 
tributions deduced by these authors are, therefore, special cases 
of the optimum distribution for a more restrictive set of condi- 
tions. 


be new, is given here. (See Fig. 1.) 


AERONAUTICAL 


This general optimum distribution, which is believed to 


| — 2x 2[ 82 V—IN) _ 
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The body, besides satisfying the fundamental conditions al. 


ready stated, is required to have volume V and areas NV at the 
nose x = 0, A at the position x = k and B at the base x = /; thys 
S S(x) dx = V, S(0) = N, S(k) = A and S(1) = B rhe opti- 


mum area distribution is 


tN + (B— N) cos™ + 2 
l 3 l 
vr /(l—x 
; ‘ 
S(x) = 
T l 16 k(l— k z — R (it — 
) (9/2? — 32 lk + 32k?) — -_ (1 — 2k) (1 — 2x) (i — 2x)24 * U ad iil : = 
{18 2 9 i? f i? 
” Ik + Ix — 2kx + 2k"? (1 — k)'/? x'/2 (1 — x)'/2) 
(k — x)? log : , 
(ik + Ix — 2kx — 2k’? (1 — k)'?2 x’? (1 — x)'?7f 
where 
(1 — 2k 1 — 2k)k'/*(1 — k)'/? 44 42(V — IN) s/t — £ 
m(A ee ee I ie : U wi i... eal = —(B- yy lF ( f) 
7 ] l 2 f 3) l ie. ip 
2 kl — k)? 
(9/2 — 32/k + 32k?) 
9 ' 
and the corresponding minimum value of the wave dray is given by 
D 4 ' j2(V — LN) it \2 ‘ (B—N)? is k(1 — k)? 
= ‘ _ — NN)? 9/2 — 3212 + 32k2)a2 
q lel l f 2 36 a ete “| 


It may be noted that a similarity law connects optimum distri- 
butions which differ only in their values of NV: 

S(x) — N =f[((V — IN), k, (A — N), 1, (B — N); x] 
D 


q 


= g{((V — IN), k, (A — N),1,(B — N)I 


Moreover, use may be made of the reversibility property that 
each optimum distribution S(x) yields another optimum dis- 
tribution by writing (J — x) for x and (/ — k) for k and inter- 
changing N and B. The area distribution illustrated is the opti- 
mum for V = 0.012/8, N = 0.006/?,k = 0.31, A = 0.014/2,and B = 
0.00972. 

Not all the parameters V, NV, k, A, J, and B may be chosen 
independently, since, for S(x) to represent a practical area dis- 
tribution, they must be such that (x) is positive for 0 < x < 1. 
However, N,/, and B may be treated as being independent: the 
value of .V is defined by the intake area, and it may be that B is 
defined by the exhaust area, while, when deducing special cases 
of the general optimum distribution, minimizing the wave drag 
alone with respect to / and B may have little significance, since 
the surface friction drag must be considered when minimizing 
the drag with respect to /, and if B is greater than the exhaust area 
the base drag must be considered when minimizing with respect 
to B. Also, k and A will always be specified together unless A is 
defined to be the maximum area C, say, when its position is deter- 
mined by the solution. Hence, the special cases which seem to 
be the most important are those for which the fixed parameters 
are V, N, C,i, and B; N,k, A,l, and B; V, N,l, and B; N, C, 
l,and B; N,l, and B. A detailed discussion of the general solu 
tion and these special cases shows that it is sometimes possible 
to achieve less wave drag than the minimum for the prescribed 
conditions by allowing V or A or both to take higher values; it 
follows that if only lower limits of V and A are given, S(x) will 


never be negative. 
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Stability of a Semiflexible Duct System* 


A. L. Guess 

Department of Aeronautical Engineering, University of Maryland 
College Park, Md 

March 24, 1954 


D™ SYSTEMS operating under high pressures and temper- 
atures sometimes find themselves being designed as column 


* This work was done while the author was employed at the Boeing Air- 


plane Company in Seattle, Wash., during the summer of 1953 
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members. Such a duct system is shown simplified in Fig. 1 
The lengths Z are circular cross-section tubes, the lengths / 
are expansion couplings, D represents the distance between 
elastic supports, and p is the pressure of the fluid flowing in the 
duct 

Neglecting pressure loss effects, Fig. 1 is essentially a ‘‘fluid 
column”’ (Ref. 1). Any axial compressive stress in the duct is 


due to the spring action of the couplings only, however, for inves 


tigating the stability of the system, the column load P = pA + 
F, must be considered, where A = internal cross-sectional area of 
duct and Fy = additional end load due to spring action of cou- 


plings at elevated temperature and pressure. The hoop stress 
caused by the internal pressure is assumed to have no effect on 
the stability of the duct system as long as it does not exceed the 
yield stress of the material used. 
Conservatively, the couplings can be considered to offer a 
flexible joint at their junctions with the ducts and, if the ducts 
L) and the couplings (/) have individually been found stable in 
the pin ended column condition for load P, the link column anal- 
ogy shown in Fig. 2 can be assumed 
considered a rigid bar pin jointed to the next member. 


each member now being 


The dashed line in Fig. 2 represents the critical deflection 
curve and the minimum energy of the system is considered in a 
manner similar to that presented in reference 2, i.e., the energy 
stored in the springs is a function of the rotation of the members 
For a system with an infinite number of members, it is found that 


(KerL/P). = 2((L/)) + 1] (Z/D)? (1) 


where 
Re 

L/D = 1 

L.i>08 


implies an infinite number of members 


= spring constant of supports necessary for stability 


Examination of Eq. (1) shows that for fixed spring constants 
K,,), the column load (P) approaches zero as L, / or D approaches 
zero. (NotethatasL~0,1S5 L/Ds This at first glance 
might not agree with physical reasoning. However, considering 
the mode of failure in Fig. 2, L or / being zero would mean that 
there could be no rotation of the members. If there is zero rota- 
tion of the members, there can be no energy stored in the springs 
and the absence of stored energy of course implies zero column 
When D is zero the springs are located at the rotation 


Hence, a rotation of the members does 


strength. 
center of the members 
not produce any energy storage and again the column strength is 
zero 

Physically though, the mode deflections that result with Z or / 
being very small might not be objectionable. If this is the case, 
a mode shape different from that shown in Fig. 2 can be assumed 
that will offer a higher buckling load. Discussion here, however, 
is limited to the mode shape shown in Fig. 2. 

The termination of the system plays a significant role if there 
are a finite number of members. For a system terminated with 
two / members, 
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Ke L/P = (KelL/P),.¢ 2 

L/2l) {(2n 1)/m)] +1 
( = { 4 

i + | 
where » = number of L members 
For a system terminated with an 1 member and an / member 

KL/P = (KaL/P).C 1) 
C, = 1/41 + (L2/D?) [1/(4m s)]] 5) 


For a system terminated with two 1 members 
KerL/P = ((K.,l./P r(D?/L? 1| ¢ (6) 
C; = 2/[(D?/L?) (27 1) +1 r 
where: + = number of / members 
The method presented above offers an approach to the semi- 


flexible duct design problem. For the ideal cases presented, the 


equations can be plotted and used quite easily by designers for 
determining a safe ratio between end load (P) and stiffness of 


supports ( K;,) 
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Further Remarks on ‘‘Interactions Between 
Wholly Laminar or Wholly Turbulent 
Boundary Layers and Shock Waves Strong 
Enough to Cause Separation’’' 


G. E. Gadd and D. W. Holder 

Aerodynamics Division, National Physical Laboratory, Teddington, 
England 

March 23, 1954 


N READERS’ FORUM FOR February, Professor Bogdonoff? criti- 
I cizes some of the experimental results, presented in reference 
1, for the upstream effect produced in a turbulent boundary layer 
by a shock wave generated by a wedge in a supersonic mainstream, 
as shown in Fig. l(a). The upstream effect is defined as the dis- 
tance between J’, the point where the plane of the shock inter- 
sects the plate, artd the point 0 downstream of which the static 
pressure begins to rise above its free-stream value po, and where 
the displacement thickness of the boundary layer is 6 Our 
experimental results (which will shortly be published in full) show 
that OJ’ /5)* is a function only of free stream Mach Number WJ 


* 


and of pg/po, where py is the maximum static pressure attained 
at the wall at the downstream end of the region of interaction 
Reynolds Number in the range covered (2 X 10% 10? based on 
the distance from the leading edge) was found to have no appre- 
ciable effect on OJ'/59* for turbulent boundary layers. We have 
not been able to obtain the report of Professor Bogdonoff and his 
colleagues,* but apparently, for high values of pp/po at My = 3, 
they found much smaller values of O/'/6)* than our experimental 
results. (In his note to Readers’ Forum, Professor Bogdonoff has 
misinterpreted reference 1: 
Fig. 1 is our experimental curve, constructed from 20-30 experi- 


the full line curve reproduced in his 


mental points, and the two points he marks as experimental are in 
fact the theoretical points calculated by the method of reference 
l We cannot accept his suggestion that the discrepancy may 
be due to our use of shock-generating wedges that did not fully 
span the tunnel, because on the contrary all those we used at 
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Numbers did fully span the tunnel. The discrepancy is all the 
more surprising because we obtained at a Mach Number of 2.5 
in a tunnel of 11 in. square working section two check results 
which were perfectly consistent with our other observations, made 
in a tunnel of 2.6-in. X 1.5-in. working section. In the 11-in 
tunnel, as in the smaller tunnel, the flat plate on which the test 
boundary layer was formed fully spanned the tunnel, as did the 
wedges used to generate the shock waves. The region of inter- 
action was about 4!/2 in. from the leading edge of the flat plate, 
so there could have been no significant interference effects from 
the side-wall boundary layers. We had originally feared that 
there might have been such interference effects in the narrower 
tunnel, but evidently this was not the case in view of the agree- 
ment with the 11-in. tunnel results. The mean curves through the 
numerous experimental points which we obtained for turbulent 
boundary layers are shown in Figs. 2(a) and (2b), in which the 
1l-in. tunnel results at My) = 2.5 are marked by crosses, and the 
results obtained by Professor Bogdonoff and his colleagues arc 
shown by the dotted curve. Fig. 2(a) shows the results for shocks 
generated externally, as in Fig. 1(a), and Fig. 2(b) shows the re 
sults for shocks generated from within the boundary layer by a 
sudden bend in the wall, as in Fig. 1(b). The upstream dis 
tances divided by 6)* are plotted against Cpr, the pressure co- 
efficient (2/yMbo?) [(pe/po) — 1]. The results of Fig. 2(b), for 
bent walls, agree well with those of Drougge.* 

In these circumstances, therefore, it is by no means evident 
that Professor Bogdonoff's results* are to be preferred to ours. 

With regard to Professor Bogdonoff’s other main point, con- 
cerning the pressure at separation, we would now agree with him 
that at 4 = 3 turbulent separation occurs at a static pressure 
ratio of only about 2.1, at high Reynolds Numbers at any rate, 
instead of 2.5 as predicted in reference 1. We had earlier thought 
that the experimental results indicated the figure 2.5 because 
dp/dx (where p is the static pressure at the wall and x the dis- 
tance along it) begins to decrease markedly downstream of 
separation when there is considerable separated flow at a pres- 
sure ratio of 2.2 to 2.5 at J, = 3. Subsequent investigations 
with Pitot tubes near the surface, however, have shown that this 
link in the pressure distribution occurs at a rather higher pressure 


than does separation. 
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Stability and Free Vibration of Nonuniform 
Beams Continuous Across Flexible Supports 


H.L. Cox and P. H. Denke 

Strength Engineers, Douglas Aircraft Company, Inc., Santa Monica 
alif 

March 30, 1954 


HERE APPEARS TO BE a necessity for treating nonuniform 
beams that are continuous over supports having deflec- 
tional and rotational rigidity. The calculus of finite differences 
affords a form of solution that may be used effectively in con- 
junction with electronic digital computers. Various boundary 
conditions may be treated without complicating the problem 
appreciably. 
The governing differential equation for a beam continuous 
across flexible supports is 


d‘ty |. d*ydI 
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dx‘ ~~ dx3 dx 


d*y d?I d*y 

, dx? 
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Fic. 1. Treatment of concentrated springs. 
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Fic. 2. Beam over flexible supports 

where 

E = modulus of elasticity 

] = variable moment of inertia of beam 

h = axial load 

k(x) = deflectional spring constant per unit length 

r(x) = rotational spring constant per unit length 

(x) = lateral load 


The sign convention used to develop Eq. (1) was the following: 


Coordinate axes positive to the right and downward 

Loads positive downward 

Shears that are positive act upward to the left of a section and 
downward to the right of a section 


Moments that are positive produce compression in the upper 
fibers 

The functions k(x) and r(x) are discontinuous when concen- 

trated spring supports exist. Actually, it is not physically pos- 


a perfectly concentrated spring; consequently, 


sible to obtain 
spring supports may be treated as shown in Fig. 1. A concen 


trated deflectional spring constant, Fig. la, may be considered 
is distributed over a finite length #. A concentrated rotational 


spring constant, Fig. 1b, may be considered as distributed in a 


triangular shape across a distance of 2h 
An example will demonstrate the general method of solution 


Consider the beam with clamped boundaries shown in Fig. 2 


Let 
a; = concentrated deflectional spring constant 
8; = concentrated rotational spring constant 
= distance between grid points 
B 
h 
. a 
k 
h 
Eq. (1) may be written in difference form by treating fictitious 
boundary points as follows.' At point (1): 
EI E 
7y1 dy2 + ys) 4 (-y 2yo + yx) (—Lo + Ie) 4+ 
: 2h4 
E p 
(—-2y + v2) (Io —- 2h 4+ kh) + —~(-2n +2) =H 
hs P h? 
ete At point (4) 
EI; E 
yo — 4ys + Gy ty; + ve) 4 -(— yo + 2s — 
2h4 
E 
25 ye) (—I3 + Is) 4+ (y 2y; + ys) (Is — 21, + Is) + 
ht 
p B 
(y 2% + ¥ — <¥ yt y 
h? 6h 
By l 
( 4 = 4; 
2h? 2h 
At point (5): 
EI E 
y ty, + 6y ty, + y (—% + ZH = 
Y ° 2h! 
E 
2¥6 + yz) (—Iy + Ie) 4+ (v3 — 2y5 + ye) (Is — 275 + Ie) + 
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etc ° 
The resulting seventeen simultaneous equations may be writ- 


ten matrically as 


ae g a p l 
D + -F+ —-G4 H + K R(L+N)|&X 
ht 2h4 hs h h 1 
Yy=Q (2 
where 
C = diagonal matrix of J values 
D = matrix of coefficients of fourth-order difference pattern 
= matrix of coefficients of second terms in Eq. (1 
G = matrix of coefficients of third terms in Eq. (1 
H = matrix of coefficients of fourth terms in Eq. (1 
AK = diagonal matrix containing a values 
R = diagonal matrix containing 8 values 
LL = matrix of coefficients of sixth term in Eq. (1 
N = matrix of coefficients of seventh term in Eq. (1 
Y = column matrix of unknown deflections 
QY = column matrix of lateral loads 
Eq. (2) may be expressed more compactly as 


(4 + = H)3 = (0 (3) 
n 


where 


E . i E P. .4 ] 
he 2h ht h hs 


The condition for buckling is 


(4 ne u) y =0 
h? 


Eq. (4) is in standard Eigenvalue form and methods of solution 
for the lowest value of p for this fundamental equation have been 
coded for various electronic digital computers 

The natural frequencies of the beam subjected to any value of 


p < Per may be obtained by writing Eq. (3) as 


(4 + pH r)s =() (5) 


where 


natural frequencies 


& 
| 


g = acceleration due to gravity 
T = diagonal matrix of beam weights 


The inertia forces of the springs are neglected in Eq. (5 It 
should be apparent that each term of the 7 matrix may be differ- 
ent. The Eigenvalues of Eq. (5) can be obtained by standard 
techniques that have been coded for electronic computers. The 
Eigenvectors also may be determined and a set of normal coor- 
dinates obtained in order to study possible conditions of forced 


vibration. 
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The Aerodynamic Forces on a Slender Low (or 


High) Wing, Circular Body, Vertical Tail 
Configuration 


Arthur E. Bryson, Jr. 

Division of Applied Science, Harvard University, Cambridge 38, 
Mass. 

March 31, 1954 


page ee THE techniques outlined in references 1 and 2, it is 
a fairly simple matter to determine the aerodynamic lift, 
pitching moment, side force, and yawing moments due to angle of 
attack, pitching angular velocity, sideslip, yawing angular ve- 
locity, and the rates of change of these quantities for the slender 
configuration shown in Fig. 1. It is more tedious but none the 
less straightforward to find the above-mentioned forces and mo- 
ments due to rolling angular velocity and its rate of change as 
well as the rolling moment due to all the previously mentioned 
motions. The only requirement is that one find the conformal 
mapping of the cross section of the configuration (see Fig. 2) 
onto a circle with no distortion at infinity. If we call the cross- 
section plane the z plane and the circle plane the ¢ plane, this 


conformal mapping is given by: 


ar 
ai 7 cr 4 1 ; w — ce® 
og 
w+te 2 oy — ce~* (1) 
C 2 ce r2 
w + =¢+ h + J -— + 
w 2 h f c 
where 
a 1 ; ; r r r ny r\2 ,A _ 
; = : sinh V2 tan 2 + Vo tan 9 + 9 tan? 9 (2) 
c/a = r/(X + sin A) (3) 
ty T 
~ - (4) 
a A sin » 
+ tan 
(h/c) + 1 (h/c) — cos X 
to T 
= . (5) 
a nN sin » 
: + tan — 
(f/e) -— 1 (f/c) + cos Xr 
r = 1/4[h + (c2/h) + f + (c2/f)) (6) 


Note that Eq. (2) 
Eq. 


and r is the radius of the circle in the ¢ plane. 
determines \ (which lies between O and 7) in terms of a/s; 
(3) determines c/a in terms of 4; Eq. (4) determines //c in terms 
of 4;/a and A; Eq. (5) determines f/c in terms of f2/a and \, and 
lastly Eq. (6) determines r in terms of h, f, and c. 

Expanding Eqs. (1) in a Laurent series which must be of the 
form: 


cert D atm 


n=0 


4c? . am . 3A cos? (A/2) 
— sin A cos* sin? - — : (8) 
3 2 2 A + sin A 


we find that: 


a, = r? — ¢c2? + 





Fic, 1. Slender low (or high) wing, circular body, vertical tail 


configuration. 
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Cross section of slender low (or high) wing, circular 
body, vertical tail configuration. 





Fic. 2. 


From reference 2, the virtual masses for motion of the cross se 
tion in the x; and x, directions are given by: 
my | F 


2rp [r? F a, (S/27)] (9 
Mize \ : 


where p = air density and S is the area of the cross section, her 


equal to ra*. Hence, substituting in the expression for a, we 


have: 
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: Fic. 3. Above: Virtual mass of the cross section for transla 
tion along its symmetry axis. Below: Virtual mass of the cross 
section for translation at right angles to its symmetry axis (4 = 
te = 0). 
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: a’ Ic? A ee 3A cos? (A/2) 
2ap|c* — — —~ sin A cos* sin*- — . 
2 3 2 4 A + sin A 
(10) 
a2 4¢2_—- e oe 
= 2rpi c* + sin A cos? sin? 
2 3 2 2 
3A cos? (A/2) 
. 4+- Dir? c?) (11 
A + sinaA 
Note that my. = O since a; is pure real. Also note that my is 


ndependent of ¢; and fz as it should be, and m2 can be written as: 


M2. = (Moe2)1,=t.=0 + 42p(r? — c?) (12) 


since yr = c when ?é, = tf = 0. Graphs of my, and m- for t; = te = 
are given in Fig. 3 as functions of body diameter-span ratio 
Utilizing mm. and mx (made dimensionless by dividing by 
oS in reference 1 and called As. and Aj) one can then find the 
forces and moments on the complete configuration as described 
inreference 1. 

From symmetry we know m;; = 0. The remaining two terms 
my and 33, which determine the roll damping and the inter- 
tions between rolling and lateral motions involve all the co- 
eflicients dy of Eq. (7), (see Ref. 2), and hence simple expressions 


for them do not seem to be readily obtainable. 
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On the Calculation of Temperature 
Distribution and Thermal Stresses in Parts of 
Aircraft Structures at Supersonic Speeds 


H. Schuh 
Division of Aeronautics, Royal Institute of Technology, Stockholm, 
Sweden 
April 1, 1954 
N TWO RECENT papers! ? analytical solutions of the temper- 
ature distribution in a skin-shear web configuration (see Fig. 
| of reference 2) have been obtained for the case of a sudden start 
{motion. However, the numerical work is rather time-consum- 
ing in any particular case. Further, in many cases the thermal 
stresses should be known for a given variation in time of aircraft 
speed and altitude, particularly for guided missiles. This re- 
quires solutions where both boundary layer temperature and 
heat transfer coefficient are functions of the time. Apparently 
such solutious are difficult to obtain by analytical means. 
Convenient and reliable methods’ are known for solving prob- 
lems of transient heat flow which use finite differences and a step 
by step procedure. The main advantage of this method is its 
flexibility and in most cases a considerable saving in time when cal- 
ulating particular cases. In the following, this method will 
e adapted to the present problem; it can be applied to cases 
where both boundary layer temperature and heat transfer co- 
ficient are arbitrary functions of the time; it can also be used 
for more complicated parts of aircraft structures than those con- 


sidered in references 1 and 2. For basic details of the numerical 


method and questions of accuracy see reference 3. 
The differential equations for a simplified configuration of skin 


nd web as treated in reference 2 read: 
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FORUM 


(7;)¢ = url Za )ez (O<x <a) (1 
(T2)¢ = xo(T2)2z2 + A(T; — T2) (a<x <b) (2) 


7,, T: are the temperatures; «x; and x2 are the thermal diffusivities 
in web and skin respectively; 4 = a/(depc), where a is the heat 
transfer coefficient, ds the thickness of the skin, p and c the density 
and specific heat of the skin; a and b — a are half the length of 
web and skin respectively; 7) is the boundary layer friction tem- 


perature. The boundary conditions are 


(7,), = 0; x = D, (T2), = 0; 


‘ 3) 
diki(T:)2 = 2deko(T2)s 


k, and ks are the thermal conductivities of web and skin respec- 


tively; d; is the thickness of the web. Initial conditions are 


T; = T2 = T; = const., 
Introducing 


1 or 2) (5) 


Pn, m9? 
Sf) 
Im = 


(g = 


TT; — T, 
Tye — Ty 
where 7, is a constant boundary layer reference temperature, 
J,, m Means a temperature at a distance x = n-Ax and at a time 
t = m:At, where Ax and Aft are the basic units in the finite differ- 
ence method. Replacing the differential equations (1) and (2) 
in the usual way by finite difference equations, the working for- 


mulas become 


Ine mei = Pi(Ant1, m + Da-1, m™) + (1 — 2pi)dn, m™ (6) 
Ono mti? = ps (Dati, mm? + Da-1, m?) 4+ 
l+s 
1— 2p2—-s s 
On, m'* (Pmt 1? + Om?) (7) 
1l+s 1l+s 
where 
x At hAt 
g , sS= , (g = 1 or 2) (8) 
(Ax)? 2 


With the help of Eqs. (6) and (7), values of 3 at a time ¢ = 
(m + 1)At are determined from those at an earlier time ¢ = mAt, 
which are known. For fulfilling boundary conditions and initial 
values see reference 3. 
of thermodynamics, the coefficients of 3 on the right-hand side 


To avoid infringement of the second law 


of Eqs. (6) and (7) must be positive. Hence 
1 — 2p, > Oand1 — 2p. —s 2 0 (9) 


determine the maximum size of p; and po. The permissible size 
of s can be found as follows: Neglect the term due to heat con- 
duction in Eq. (2), solve the equation in a closed form and com- 
pare it with the corresponding solution of the finite difference 
method. If # and 7; are constant and ¢ is the permissible error 


in 7, — 7, at a certain time ¢, the condition for s is 


Dm” ht/2s 
1 — (; ; eM! < |e 
+ § 


As to the choice of the size of Ax see reference 3. 
value of Ax is chosen, the number of steps in the computing 
scheme will be the smaller the bigger p is; hence it is convenient 
to choose At as big as possible, subject to the inequalities (9) and 
Assuming the same material for skin and 


(10) 


If a certain 


(10) being fulfilled. 
web, Af is thus obtained from (9) as 


2« hy, i 
At = + (11) 
(Ax)? 2 
where h,, is the maximum value of A, if # is variable. For the 


skin-shear web configuration of reference 2, three particular 
cases have been treated, whose steady state is the same; the vari- 
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ation of other characteristic parameters with time can be seen 
from Fig. 1. The cases are: (i) sudden start of motion; (ii) 
linear increase of boundary layer temperature, with heat transfer 
coefficient assumed constant; (iii) linear increase of Mach Num 
ber in horizontal flight (i.e., constant acceleration) with vari- 
ation of boundary layer temperature as well as heat transfer 
coefficient under the (a) turbulent 
boundary layer on the skin, hence h ~ .\/°* and the temperature 
recovery factor r = a /3, Prandtl Number = 0.7 
(air), (b) Amar equal to the value of / used in (i) and (ii). Cases 
(i) and (ii) have been chosen for comparison with analytical 


following assumptions: 


where o¢ = 


solutions;' ? a case similar to (iii) has apparently not been treated 
before. For all three cases Ax = 1.17 in. and other values have 
been assumed as in reference 2. At follows from Eq. (11). The 
temperature distribution for (i) obtained by the present method 
is compared in Fig. 2 with those of reference 2. Agreement is sur- 
prisingly good, considering that half the web is divided into four 
intervals only and that the computation requires only 2 hours 
(for (iii) about 3 to 4 hours). 

0 (middle of the web), are shown in Fig. 3. 
3 per cent for 


The greatest stresses, which occur 
atx = For the 
maximum stress the deviation is not more than 
(i) and somewhat bigger for (ii), where the analytical solution is 
thought to be less accurate, particularly at small times; besides 


the stresses of (ii) should have a horizontal tangent at ¢ = 0. 


— Analytical method + xy 
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